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Abstract 


In this thesis we start by reviewing theoretical aspects of micromagnetism. Since many tech- 
nological applications only depend on the behavior of magnetization of ferromagnets with the 
applied magnetic field at a macroscopic scale, there is no need to use a highly detailed theory 
like quantum mechanics. Micromagnetics is the framework that captures at a mesoscopic 
level the essential dynamical behavior of a magnetization field. It describes the combination 
of a very fast processional motion and a slower damping toward the magnetic field. 

The two central relations are the Landau-Lifshitz (LL) and Landau-Lifshit-Gilbert (LLG) 
equation’s. We show how to derive the first assuming two main experimental and theoretical 
observations: (1) the local magnetization norm is conserved; (2) the equilibrium state that 
both the magnetic field and magnetization aligned. We then analyze the dynamics implied by 
four magnetic held contributions: the applied held; the anisotropy held which has a similar 
behavior to an applied held along the lattice axis; the stray held which depends on all other 
magnetic moments; the exchange held which is the most relevant term, it tends to smooth 
out the magnetization direction. 

We then introduce the LLG equation by representing the damping term by Rayleigh 
Dissipation. It is an implicit equation of the magnetization. Our goal is to develop a Python 
code to integrate this equation. We start then by combining it with the Finite Element 
Method to discretize space and the Implicit Midpoint Rule to discretize time. To avoid 
meshing the surroundings of the system that was required for introducing the asymptotic 
boundary conditions for the calculation of stray held potential, we use the Boundary Element 
Method and a new potential to restrict these calculations to the system. Using the Newton 
Raphson Method we obtain a linear system of equations that is solved at each time step 
yielding the evolution of the magnetization. 

Setting our code to solve a standard problem for permalloy block 120 x 120 x 10 nm 3 we 
compare our results of magnetization evolution with those of the OOMMF micromagnetic 
simulator to validate our code. The results of our code compare very well with those of the 
OOMMF simulation, specifically the time evolution of y component of magnetization, its 
Discrete Fourier Transform as well as the spatial distribution of the amplitudes of the Fourier 
coefficients for two distinct resonance frequencies. 


Keywords: Micromagnetism; Landau-Lifshitz; Landau-Lifshitz-Gilbert; Rayleigh Dissipa- 
tion; Implicit Midpoint Rule; Finite Elements Method; Newton-Raphson Method; Boundary 
Element Method; Discrete Fourier Transform; Python 
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Resumo 


Ao nivel quantico, os materials ferromagneticos sao caracterizados por magnetismo resultante 
do spin dos eletroes cuja descrigao completa requer a aplicagao da Mecanica Quantica. Estes 
materials quando TcT c sao caracterizado por dominios i.e. regioes de momentos magneticos 
com a mesma direcgao e sentido, que sao o resultado de fortes interagoes de troca. A orientagao 
da magnetizagao local varia a escala destes dominios quando o sistema esta ern equilibrio. 
Quando um carnpo magnetico externo e aplicado sobre o sistema, a magnetizagao de cada 
dominio tende a alinhar com este campo. Se o carnpo for suficientemente forte passamos de 
multiplos dominios, a um unico monodominio. Quando este campo e desligado, estes voltarn a 
formar-se mas com orientagbes diferentes das iniciais. Este comportamento da magnetizagao 
total do ferromagnete em fungao do campo aplicado designa-se por curva de histerese e a 
sua taxa de variagao quando o campo e zero determina a estabilidade da magnetizagao total. 
E o formato desta curva que permite aplicagoes tecnologicas de gravagao e leitura de dados 
em discos rigidos, construgao de eletroimanes ou libertagao de calor para destruir celulas 
cancer osas. 

Com o objetivo de aplicar estes materials e preciso tragar a curva de histerese. Observa-se 
que e uma relagao entre variaveis na escala macroscopica. Conclui-se assim que uma descrigao 
quantica do material seria desnecessariamente rigorosa. O que e necessario e informagao a uma 
escala mesoscopica que capture as caracteristicas observadas quanticamente e que reproduza 
o comportamento magnetico a uma escala macroscopica. A area de trabalho que estabelece 
esta ligagao e a de Micromagnetismo. Ao inves de partir dos momentos magneticos associados 
a electroes, divide-se a amostra em subvolumes que sao grandes quando comparados com 
a constante de rede do ferromagnete, mas rnuito menores que o tamanho da amostra. A 
informagao de corno a magnetizagao varia no interior destes subvolumes e desprezada. Em 
relagao a dinamica da magnetizagao duas equagoes foram desenvolvidas uma por Landau e 
Lifshitz, e outra por Landau, Lifshitz e Gilbert. Estas equagoes descrevem a combinagao 
de dois movimentos do vector magnetizagao e dependent do campo magnetico local. Um 
movimento rapido de rotagao da magnetizagao em torno da direcgao do campo magnetico e 
outro movimento mais lento de amortecimento que descreve a perda de energia do sistema. 
A evolugao temporal deixa de ocorrer quando a magnetizagao esta alinhada com o campo 
magnetico local, obtendo-se finalmente uma configuragao mesoscopica que esta associada a 
um ponto da curva de histerese. 

Nesta dissertagao irernos rever primeiro como obter a equagao de Landau-Lifshits (LL) a 
partir de duas suposigbes: (1) a conservagao da norma do vector de magnetizagao local; (2) 
o estado de equilibrio e atingido quando magnetizagao e campo estao alinhados. A seguir 
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apelando a mecanica classica, Gilbert propos representar o termo de amortecimento via dis- 
sipagao de Rayleigh obtendo-se assim a equagao de Landau-Lifshits-Gilbert (LLG). 

Nestas duas equagbes o campo magnetico local apresenta quatro contributes: campo apli- 
cado; anisotropia; campo desmagnetizante; campo de troca. No equilibrio final observa-se o 
alinhamento entre magnetizagao e campo magnetico. 

Ambas as duas hipoteses acirna sao validas para todos os contributes magneticos, a 
diferenga fundamental so esta na fonte do campo. O campo aplicado e determinado ex- 
ternamente. O termo de anisotropia e representado por um campo ao longo da estrutura 
cristalina do material. O campo desmagnetizante resulta da soma dos campos de todos os 
outros momentos magneticos. O termo de troca ferromagnetico tende a orientar momentos 
magneticos vizinhos no mesrno sentido. 

Comegamos com a equagao de LLG, uma equagao implfcita na magnetizagao, e o campo 
desmagnetizante determinado pela equagao de Poisson. Sobre o sistema de equagbes LLG e 
equagao de Poisson aplicamos o “Implicit Mid-Point Rule” , discretizando as derivadas tempo- 
rals sob um certo incremento de tempo e substituindo a magnetizagao e o potencial por rnedias 
entre estes dois instantes de tempo. Obtem-se assim duas equagbes que os relacionam e que 
correspondem a um sistema implicito sobre o proximo estado ern fungao da magnetizagao e 
do potencial no estado anterior ou inicial. 

Dada a complexidade destas equagbes e o facto de o incremento temporal ser pequeno, a 
magnetizagao e o potencial nao variam muito neste intervalo, donde deixa de ser necessario 
informagao sobre variagoes de ordern superior a primeira. Expandindo os residuos destas 
duas equagbes ate primeira variagao do campo e do potencial, obtem-se como coeficientes os 
Jacobianos dos residuais. Este e o metodo de Newton- Raphson e e muito usado para resolver 
equagbes com escalas de tempo muito diferentes, como observado nas equagbes LL e LLG. 
Com este rnetodo define-se o residuo associado a equagao LLG e a equagao de Poisson e 
procura-se a solugao do sistema de equagbes correspondente. 

Para implementar computacionalmente a resolugao destas equagbes que sao contfnuas 
no espago e preciso discretiza-las utilizando-se o Metodo de Elementos Finitos. O sistema 
e a sua vizinhanga sao divididos nurna rnalha de tetraedros e a cada um dos vertices de 
cada tetraedro associa-se um vector com as tres componentes de magnetizagao e o potencial 
magnetico. Supondo que sao fungoes lineares e possivel representa-las como combinagao 
linear de fungoes de base tambem lineares e que dependem da localizagao dos vertices de cada 
tetraedro. Desta forma o sistema de equagbes dos residuos pode ser simplificado, mas devido a 
presenga de segundas derivadas espaciais nos terrnos de troca e de “Poisson” e essencial reduzir 
a ordern destas equagbes. Para tal usa-se a forma fraca do sistema de residuos expandidos em 
primeira ordem. O problema inicialmente continuo passa a um sistema de equagbes lineares 
em que a incognita e um vector cujas componentes sao todos os incrementos de magnetizagao 
e potencial em todos os vertices da rnalha e com o qual e possivel actualizar a configuragao 
em cada instante. 

Para evitar ter de usar uma rnalha para o reservatorio (espago no exterior a amostra) ire- 
mos introduzir um novo potencial que satisfaz uma nova equagao de Poisson que tern condigoes 
de fronteira de Neumann na fronteira do sistema. Com o Metodo Elementos de Fronteira e 
possivel rnapear a solugao deste problema na fronteira para o potencial do campo desmagne- 
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tizante tambem na fronteira. Desta forma conseguimos substituir as condigoes de fronteira 
assimptoticas que originalmente requeriam a presenga do reservatorio, por um problema de 
Poisson com condigoes de fronteira de Dirichlet. O custo desta reformulagao do calculo do po- 
tencial do carnpo desmagnetizante e o aumento do sistema de equagoes a resolver e o calculo 
da rnatriz que rnapeia estes potenciais, que e densa e computacionalmente exigente de ser 
obtida. Contudo as suas entradas sao constantes donde este bloco do sistema de equagoes so 
precisa de ser calculado uma vez. 

Cada um dos outros termos do campo magnetico tambem ira contribuir com o seu Jaco- 
biano. Apos introdugao de uma quadratura nodal para aproximar os integrals ai presentes 
obtem-se a forma final para ser implementada computacionalmente. 

Com estas expressoes e os dois problemas de Poisson e possfvel calcular as entradas da 
rnatriz. Dado que tal tern de ser feito a cada incremento temporal o processo de calculo 
e construgao da rnatriz tern de ser rapido. Daf que tenhamos desenvolvido os algoritmos 
em Python para cada terrno de LLG usando operagoes vectoriais. O sistema de equagoes 
pode agora ser resolvido o que permite actualizar um estado corrente de magnetizagao e dois 
potenciais ao longo de todo o sistema magnetico para um estado seguinte. Constroi-se a 
historia do sistema magnetico, donde calculos da sua curva de histerese sao possiveis e assim, 
finalmente desenvolver as suas aplicagoes. 

Na parte final deste dissertagao para validar o nosso codigo implementamos a resolugao de 
um problema standard sobre um bloco de permalloy de 120 x 120 x 10 nrn 3 . Obtivemos resulta- 
dos para a evolugao da componente y da magnetizagao e da Transformada de Fourier Discreta, 
assim como a distribuigao espacial das amplitudes e fases associados a duas frequencias de 
ressonancia distintas. Estes calculos sao comparados com os resultados publicados para o 
mesrno problema, concluindo-se que estao em boa concordancia. 


Palavras-Chave: Micromagnetismo; Landau-Lifshitz; Landau-Lifshitz-Gilbert; Dissipagao 
de Rayleigh; Metodo Elementos Finitos; Metodo Newton-Raphson; Metodo Elementos de 
Fronteira; Transformada de Fourier Discreta; Python 
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Chapter 1 


Introduction 


At the fundamental structure of materials we observe the orbital angular momentum of elec- 
trons and spin angular momentum of electrons [1]. Orbital magnetism is associated with 
diamagnets and superconductors , while electron spin is associated with paramagnets and fer- 
ronragnets [2], In classical electromagnetism all of these are uniquely represented by a field J, 
that is interpreted at every spatial point x as indicating how much charge crosses a given area 
element per unit time [3]. The way these charged currents affect themselves is represented 
by a magnetic field, H. The framework that describes this interaction is in such a way that 
other representations of these currents are more adequate. Instead of charged currents, it is 
possible to define the magnetization field, M, that also takes into account the geometry of 
the currents. 

The way charged particles affect each other is reformulated as describing how the magne- 
tization field in one region affects other regions. The dynamics that comes as a result of this 
interactions we call the process of magnetization of materials. 

All materials exhibit magnetization at some scale. In this thesis we will focus on a subset 
of these materials known as ferronragnets. These respond to externally applied magnetic fields 
in more complicated manner than the simple relation M = \'H ap , and have a relation not 
represented by a closed function but by a non-singular and multivalued relation represented 
by a plot and designated curve of hysteresis , [3, p. 421, 443]. 



Figure 1.1: Hysteresis curve for a ferromagnetic material represents how the total magnetization M of the 
sample behaves as we change the applied field, [3]. 
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It is a macroscopic property that can be experimentally obtained with great precision 
using a Superconducting Quantum Interference Device [4] . The macroscopic magnetization of 
the entire sample is given as a function of applied field. The causes of such behaviour can 
be understood using a mesoscopic approach. Introducing the existence of magnetic domains 
in the sample, that are subparts of the system where magnetization direction and norm are 
independent of x, the variation occurs only between domains. The existence of magnetic 
domains can be seen in the following way: Exchange interaction drive the entire sample into a 
single domain, but if we introduce interactions like stray field, this single domain is broken in 
smaller ones. The corresponding field can be treated analogously to an electric field, defining 
fictitious magnetic charges, with volume density given by p* = —V • M(x) and a surface 
density given by a* = M(x) • h where M(x) is the magnetization at x. In this way the 
monodomain in figure 1.2 second picture would have p* = 0 in its interior but the top surface 
would have ’’positive” charge while the bottom have the opposite. 



Figure 1.2: Two domain configurations. The first one has a lower energy [3]. 

Apelling to arguments from statistical physics, the configuration in figure 1.2 first picture 
would have a lower energy since the charges cancel at domain boundaries and M is orthog- 
onal to the surface of the system, hence energy associated with stray field is zero and this 
configuration is the one observed. However in real materials, it is observed that the existence 
of these domains also depends of the presence of strong anisotropy effects in the material [5] , 
i.e., preferable directions of alignment within the material. 

In figure 1.1 we see at the beginning of the path, when all domains are disaligned, |M| = 0, 
corresponding to first picture of figure 1.2. This configuration minimizes the sum of exchange 
energy and stray field energy. 

However if we introduce an applied field, the magnetization will tend to align with this field 
H ap . The domains in line with this field grow progressively at the expense of the disaligned 
ones as they rotate toward the magnetic field. After a given intensity, the sample is only one 
domain, second picture in figure 1 . 2 . 

Now, if the field is reduced to zero, the magnetization will not retrace through the previous 
path, instead a new state is reached where |M r | 7 ^ 0, a remanescent magnetization prevails 
like in permanent magnets. Domain emerge again but not the same as in the beginning of this 
path. The new magnetization configuration retain some of the properties that characterize 
recent states, in this case the state where all magnetization was aligned with the applied field. 
Continuing applying the magnetic field in opposite direction, the magnetization will reach 
zero at — H c j and finally if we continue increasing the magnetic field in the opposite sense, 
the entire sample will align in the opposite direction. Inverting the process we can close the 
cycle. 
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Notice during the entire process the magnetization norm, |M(x)|, inside each domain is 
constant. 

The behavior of domains shown by the hysteresis loop, allows very important applications. 

• Consider then the following simple idea. With a small set of ferromagnetic samples. Each 

one will have their own hysteresis curves and under an applied external field each one will 
behave according to it. If we magnetize one sample in the y direction, after turning off 
the field, the monodonrain disappears and several domain appear but in average they will 
point in the y direction. This configuration saved information in the most elementary 
form, we can call it 1. If we assume that magnetizing in the opposite direction means 0, 
then we can easily see that with a set of these materials we can encode information in 
binary format. The information can later be retrieved. Since magnetic moments point 
on average in one direction previously decided, they generate a magnetic held in that 
direction which can be measured. For example with materials whose resistance varies 
with the direction of the held, the held direction can be converted into one of current 
intensity within the measuring device. 

This is the main idea behind hard drive recording and reading, the direction of the mag- 
netic moments on the hard drive encode the information we introduce. The guarantee 
that this information is saved through time can also be seen in the hysteresis curve. If 
we have an high H c i , the curve will be very wide, so the derivative of magnetization 
with the applied held at H = 0 is close to zero, for the other small helds present within 
the hard drive and even external ones the magnetization will not change much, and the 
information is preserved. 

• Ferromagnetic materials are also used in boosting the intensity of magnetic helds. Suppose a 

coil of electrical wire. Under a given current, a magnetic held is generated, if at the center 
of this coil this material is introduced, then it will be magnetized, the magnetic held 
produced by this sample will add to the held produced by the coils. The electromagnet 
is turned on to attract and must reach — H c i to the electromagnet be off. 

Since the energy spent is given by 

Sw = n o J H ext ■ 5Md 3 x (1.1) 

V 

we can see that if ff c i is close to zero [3], the energy cost is lowered. A thinner hysteresis 
curve is cheaper! 

• If to the work necessary to magnetize a sample to <5M, we subtract the variation of exchange 

energy, we obtain the heat produced by the material [6, p. 180] 

<5Q = do {H ex t + n w M)5M (1.2) 

Aligning the domains imply heat production. This can be used with medical purposes 
to eliminate tumors. Developing particles with receptors for these cancerous cells, after 
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they attach we apply an external magnetic field. The heating produced by the particles 
is transferred to the surrounding cells destroying them [7]. 

These examples show us not only that initial and final states of equilibrium magnetization 
are important, like we also saw in the hysteresis loop, but the dynamics of relaxation is 
important as well if we want to understand how to design materials with new hysteresis 
shapes. 

Still within a mesoscopic framework we can build a theory that explain the dynamics 
of magnetization. Proposed by Landau and Lifshitz in his paper [8] and later reformulated 
by Gilbert in [9] these phenomenological equations constitute the foundations of Micromag- 
netisrn. 

In this thesis we first introduce the theoretical background for Landau-Lifshitz equation 
in Chapter 2 and for the remaining Chapters we intent to develop the code to integrate 
the Landau-Lifshitz-Gilbert equation by applying the Finite Element Method and Newton- 
Raphson Method giving us the magnetization dynamics. 
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Chapter 2 


Landau-Lifshitz and 
Landau-Lifshitz- Gilbert Equations 


Micromagnetism is the theory that has the goal of using classical electromagnetism to model 
the behaviour of magnetic materials. As the name suggests, it describes the material at the 
micrometer scale. 

On the atomic scale, magnetization is a result of charged currents associated to spin and 
atomic orbital moments whose quantum description is highly complicated. However approx- 
imations can be made if the final output of the theory does not require all the microscopic 
details. For ferromagnetic materials it is known that at T <C T c , the system’s magnetic 
structure is characterized by domains of aligned magnetic moments due to strong exchange 
quantum interactions that with crystaline anisotropy have the tendency to be oriented along 
some lattice directions of the sample [2]. 

Under these conditions a detailed quantum description is unnecessary, it is then open a 
strategy of coarse graining that has the goal of simplifying the description of the magnetic 
moments and the dynamic laws they obey [11]. 

In what follows the dynamics of an out of equilibrium system of magnetic moments sub- 
jected to a given magnetic held will be derived from two fundamental principles giving us 
LL equation, then from energy considerations, the magnetic held contributions are obtained 
through variational calculus, the final equilibrium point condition is derived. With the LL 
equation each of the magnetic held contribution has dynamical implications that will be anal- 
ysed. Finally a second formulation of LL equation is introduced, Landau-Lifshitz-Gilbert 
equation, which describes a dynamics similar to the hrst and that we will use in subsequent 
Chapters. 


2.1 Landau-Lifshitz equation 

Each volume element is characterized by a current j m which is described by quantum me- 
chanics and is due to the motion of charges [3, p. 408-410]. To simplify we define M(x)H to 
be the average of magnetic moments in the vicinity of x in a cubic cell of volume, H, of side 
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much bigger than the lattice spacing but much smaller than the sample size, it is given by 

M(x) = = 2^ / S X ' iM C 2 - 1 ) 

n 

where m is the total magnetic moment associated to volume f i. 

Once the magnetization is known, the details of how it is generated are irrelevant. Despite 
having an explicit relation, its computation will not be needed. It is a continuous function that 
will not have variations on scales of the order of lattice spacing, the details of how variations 
occur within the sample subdivision don’t belong to micromagnetism, instead the theory will 
focus in determining how variations in M on the scale larger than this sample subdivision 
evolve in time when the system is out of equilibrium. 

Following the coarse graining strategy (CG) used above for magnetization, its fundamental 
dynamical processes can also be recast in this new scale. Despite the fact that j m comes from 
spin and orbital magnetic moments [3, p. 409], when we use m or M defined in 2.1 both this 
sources are irrelevant, therefore to understand the dynamics of m we can use a simple example 
from classical electromagnetism. Suppose a current distribution smaller than a lattice cell and 
that satisfies V • j = 0, for example a circular loop of current. The following identity then 
follows 

V • (x k xa) = xij k + x k ji (2.2) 

where k and l indices indicate the components of the position vector x = (x,y,z) T . 

Using the divergence theorem we find 

J d 3 xj k xi = - J d 3 xj k x k (2.3) 


from which we can derive 


J d 3 xj k xi = ^ J d 3 x(j k xi - jix k ) (2.4) 

If the current j is subjected to a magnetic field B = /io(M + H) then the torque is given by 


N = / d 3 x x x (j x B ) 


(2.5) 


We would like to reformulate it in terms of magnetic moment defined in (2.1). Since x and j 
are already present we only need to join them. Using the identity 


a x (b x c) + b x (c x a) + c x (a x b) = 0 


we get 


N = — j d 3 x j x (B x x) - J d 3 . x B x (x x j) 

Introducing now a x (b x c) = b(a • c) — c(a • b) yields 

N = — [ d 3 x B(j • x) + [ d 3 sx(j • B) — [ ci 3 xx( B-j)+ [ d 3 xj(B-x) 


(2.6) 


(2.7) 


(2.8) 
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The first integral will be zero for a circular loop of current as j will be perpendicular to x. 
For arbitrary shape current distributions, since we assumed they are small, we can consider 
B constant and factor it. The first integral can then be shown to be also be zero using (2.2) 
and the divergence theorem 

J d 3 xj x x + j y y + j z z = ^ J dS ■ j(x 2 + y 2 + z 2 ) = 0 (2.9) 

where S encloses the current and therefore j = 0. The torque is reduced to the last integral 1 
which is reformulated by components in terms of the magnetic moment using (2.4) and where 
factoring B yields 

N k = B t J d 3 xj k xi = Bi^ J d 3 x ( j k x t - j t x k ) = e kH Bimi (2.10) 

These components define now the torque as 

N = m x B (2-11) 


We can relate the torque with magnetic moment by first relating it with the angular momen- 
tum. Suppose the current j = q k v k 5(x. — x k ), then substituting in (2.1) gives us 

m = 1 J>(x t X v t ) = £ ^ L ‘ = ^ L (2.12) 


where the angular momentum of the particle k is given by L k = x k x (r nv k ) and the last step 
assumes all k particles are equal and the total angular momentum is L = h k . This simple 
relation allows the reformulation of (2.11) in terms of the magnetic moment of the current 


dm 

dt 


7 m x B 


(2.13) 


where we define the gyromagnetic ratio 2 as 7 = ^ . If initially the current magnetic moment 
is not parallel with the magnetic field and q > 0 then it will precess around it in a clockwise 
way and will never align with it. Observing now that 4 and xB operations are linear, our 
CG strategy can again be applied. Since a current distribution is described by this equation, 
so each one of the many current in a volume 12 will have dynamics described by a similar 
equation, whose net effect is 



3 


dnij 

dt 


1 

n 


7m x b 

3 


d h E m j 

3 

dt 


= 7 



dM 

dt 


7M x B 


x B 


(2.14) 


1 Observe that (2.4) can be used to show that f rf 3 ix( B ■ j) = m x B. 

7f q is the electron charge, q = \e\, then 7 — > j e , the electron gyromagnetic ratio. 
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which can be rewritten in terms of H using B = p o(M + H) giving us 


= 7 /^oM x H (2.15) 

At the scale of each cell the magnetization M evolves in time by precessing around the local 
magnetic field as well. However it is a experimental fact that eventually the magnetization 
will align with the field, and also that for temperatures below T c its norm is conserved [3, 11]. 
Rather than derive a damping term from fundamental quantum processes, with no more 
knowledge than these two principles we can look for the simplest dynamical expression that 
would be compatible with them [11], and then test its consequences against experiment. 
Suppose then a first order dynamical equation that depends on the magnetization and the 
held, H: 

^=V(M,H) (2.16) 

Given that the norm is constant, there is no variation along the magnetization vector, so 
it must be an important direction, hence we construct an orthogonal reference frame that 
includes it: {M,M x H.M x (M x H)}. By representing V in it, the equation becomes 


— — = aoM + aiM x H + a^M x (M x H) (2-17) 

Since the norm of M does not change in time, then a 0 is zero. The second principle states that 
an equilibrium magnetization configuration points along the magnetic held, that is M e(J x H = 
0, but that is readily satished by the equation because then ^ = 0. 

What is left is the determination of ai and 02 , which must include po as we are expressing 
the dynamics equation in terms of H instead of B. The sign of the hrst must negative if the 
charge is q = —e in (2.15), the second must also be negative by noting that then this second 
term updates the magnetization towards the magnetic held. From (2.15) the constant ai will 
be the electron gyromagnetic ratio times the vacuum permeability, 7 l = 7 e/Uo, and it controls 
the rate of precession of magnetization, it has units The constant <32 is decomposed as 
[11, p. 26], where a is a small dimensionless damping constant usually between 10~ 4 — 10~ 2 . 

The equation obtained is the Landau-Lifshitz equation 


dM 

dt 


— 7lM x H 


— M x (M x H) 
M s y J 


(2.18) 


It describes the continuous dynamics of a magnetization function subjected to a magnetic 
held, characterized by precession about H and damping that represents the thermal contact 
with the reservoir. Notice however that nothing was assumed about the source of such held. 
Several contributions will be taken into account: an externally applied held H ap generated by 
the reservoir; the influence of lattice structure directions, H/yv; the demagnetizing held H m\ 
and the microscopic representation of quantum exchange The total magnetic held in 

LL equation (2.18) is 


H tot — H ap + Hyuv + Hj\/ + H ex 


(2.19) 


In the next section each contribution to the dynamics and energy will be analysed individually. 



2.2 Energy and Dynamics of Magnetization 


The energy associated with the system under a given field contribution will be given. For both 
the external field and stray field the energy is readily obtained from the fields, for anisotropy 
and exchange, the energy function form is first introduced and from it, through variational 
calculus, its field representation is obtained. The implication of each magnetic field on the 
dynamics will be studied on next sections. But first it will be useful to simplify all the 
equations by identifying the natural scales for each quantity, it is possible then to put both 
the LL and energies in adimensional form. 


Fundamental scales in Micromagnetism 

The dynamical equation (2.18) contains parameters with units, with these and its combi- 
nations we can construct the fundamental scales for this particular problem. Re-expressing 
all quantities as adimensional will isolate its mathematical form, which has no units, and is 
therefore more simple. 

The magnetization, M, has norm M s and units so it can be written as M = M s m, 
where the new field 3 m point in the same direction as the original, but now its norm is 1. 
The magnetic field whose norm can change has the same units as magnetization and will be 
rescaled by M s as well, H = M s h. The spatial characteristic length has to be build by a 
combination of variables. Knowing that exchange interaction is characterized by a parameter 
A of units Jm _1 and that yoM 2 has units of Jm~ 3 , then defining Iex = \J fl ^M 2 have 
units of m, where the /M) is the vacuum permeability. Therefore we can measure lengths as 
x = ul ex • For a generic function f(x) spatial derivatives will be given by 


df(x) _ df{ul E x ) _ df(x) dx 
du du dx du 

Which can be rewritten as 

df(x ) _ 1 df(ul E x ) 
dx l ex du 

Gradients and Laplacians are composed by derivatives, then 


( 2 . 20 ) 

( 2 . 21 ) 


Vu/Og y, z) = lBiV x /(r, y, z) (2.22) 

Vu/(x, y, z ) = l 2 EX Vlf(x, y, z) (2.23) 

From the units of /to AT 2 , introducing now the volume l 3 E xi energy scale is given by E = 
e/io M 2 lg X , where e is dimensionless. Any energy term can be put in this form. 

As we shall see in Chapter 3, the demagnetizing field has a potential formulation that 
satisfies a Poisson’s equation V 2 </> = V • M which must be adimensional as well. Defining 
the unit of potential, <F, we have cj) = (f> ad & . Substituting in Poisson’s equation together with 
(2.22) and (2.23) we have 

JLv 2 ^ = • m (2.24) 

'-EX EX 

which suggests that the fundamental unit of potential is 4> = M s /Iex- 
3 This m is not the magnetic moment we defined in (2.1). m is now defined as 
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Finally the fundamental unit of time is obtained by combining the precession constant 7 l 
with the magnetization norm, M s , giving us ( 7 lM s ) -1 = ('y e /-io M s ) _1 with units of s. 

With these scales, substituting the (2.21) version for time unit into the LL equation (2.18), 
and factoring the fundamental units out of magnetization and magnetic fields the adimensional 
LL is 

( 1 1X1 

- — = — m x h — am x (m x h) (2.25) 

at 


Not only it is cleaner, it brings attention to the important fact to the magnitude difference 
that separates the precession and damping terms. While the first has coefficient 1, the second 
is determined by a, which is 2 to 4 orders of magnitude smaller, the rotation motion is 
much faster than the damping. So huge is the difference between this two behaviour that 
differential equations that exhibit this kind of difference have their own name, stiff equations. 
Each contribution for h is now analysed where we assume h, m and x are dimensionless. 


Contribution 1: The field generated by the reservoir 

If a magnetic material is subjected to an applied field h ap then the energy will be given by 

e ap = - J h ap ■ m d 3 x (2.26) 

v 

The energetic formulation allows to anticipate which field configurations m are the most 
benefited. A microstate of the universe specifies the state of the sample, m, and the reservoir. 
The energy of the universe is constant hence for a given m there will be several reservoir 
configurations compatible with the total energy e un i verse = e s i s t + e reserv . Hence for a given 
microstate for the system with energy, e s i s t , the number of microstates of universe is pro- 
portional to the number of states of the reservoir, £l'{e S i S t). It is known this quantity is 
proportional to e _/3fW . It is expected that if e s i s t = e ap given above, then the more aligned is 
the magnetization with the applied field , the lower is e ap and the higher will be the number 
of reservoir states, £l'(e s i s t). This system configuration is therefore more probable when the 
system is in equilibrium. It is in this sence that the energy associated with a given field is said 
to contribute or penalise certain magnetization configurations of the system [13]. Despite the 
fact this argument is for externally applied fields, the same will work if the source of the field 
is other. 

The dynamics described by the applied field can be understood discretizing time in LL 
equation giving us 


rrL = m — m x h At — am x (m x h) At (2.27) 

where m' = m(f + At) and m = m(t) and h = h ap . From the first term, at each time step it 
always updates the current state m with an increment m X h At that is orthogonal to m. If 
At — > 0, this dynamics is represented by a circular trajectory around h ap . 

The second term will add an increment that in turn is orthogonal to the first term, meaning 
it evolves m into pointing along the applied magnetic field. Observe that as this alignment 
occurs |m x h| gets smaller since |m| = 1 and the applied field is constant. Therefore the 
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precession motion becomes gradually slower until it ceases when magnetization and the field 
are parallel. 

The dynamics associated with h. ap describes for many steps, the process that occurs at 
one step when the field on the sample is not the applied field, h ap , but a more complicated 
field that at each step is changed, it can change because it will depend on the current system 
configuration, m(x, t). These are the next contributors to be analysed. 


Contribution 2: Uniaxial anisotropy 


The lattice structure of a magnetic material benefits some spatial directions for the magnetic 
moments to align. A rigorous description would require quantum mechanics, but a simple 
formulation can be made the following way. Suppose that m belongs to a magnetic domain 
whose lattice has e y as a favoured orientation. To describe this phenomenon it is introduced 
a magnetic field in that direction given that as in, h ap , a field in that direction will make 
m evolve until it coincides with e y . However, unlike an external applied field, a favoured 
direction only requires that m evolves until it coincides with the direction of e y , not its sense. 
To describe this, it is introduced on this e y held a coefficient that depends on magnetization, 
n ad (e z ■ m), where the adimensional constant for uniaxial anisotropy, K ad , is given by fl 2 £p 
[11, p. 33]. When positive, m evolves until superposition with e y , just like an applied held, 
h ap , in this direction, but if (e z ■ m) is negative, m evolves until it coincides with — e z . 

For a generic favoured direction, this held is then 

h A N = e AN K ad (e AN ■ m) (2.28) 


The energy formulation of such held is given by 

e AN = J /Av(m) d 3 x (2.29) 

V 

K ad 

Ian( m) = -^-(1 - (rn • e AN ) 2 ) (2.30) 

Observe that the energy is minimum when m coincides with e AA r, that is (m • e^w) 2 = 1- 
These are the most benehted configurations. To show the definition of f A n leads to (2.28) we 
can recover the anisotropic held from the energy. Through variational calculus, suppose we 
have the hrst order variation in m is given by 


ad r ^ ad r 

5e AN = — I (1 - ((m + (5m) • e AN ) 2 )d 3 x — I (1 - (m • e AN ) 2 )d 3 x (2.31) 

K ad f f 

= — / (dm ■ e AN ) 2 d 3 x - K ad / (m • e AN )(Sm ■ e AN )d 3 x 


The squared variation makes the hrst integral much smaller than the second, reformulating 
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the latter we have 


Se A N = — n ad J [(m • e A N)e A N] • Smd 3 x 
= -n ad [ h AN dm d 3 x 


(2.32) 


where we identify the anisotropy field (2.28). 


Contribution 3: The magnetic field of magnetization 


The demagnetizing field, h m, can have a potential formulation satisfying a Poisson’s equation 
under appropriate boundary conditions which will be the subject of next chapters. For now, 
its explicit solution is given by 


h m = 


1 

47T 



X — X 


X — X 


/|3 



dS' ct*(x / ) 


X — X 


X — X' 


./|3 


(2.33) 


where p* = — V • m and a* = m • h are interpreted respectively as the density of fictitious 
magnetic charges within the system and on its boundary in an analogy with the electric field 
from polarized dielectric materials [3, p. 416]. 

The sum of energies associated to each magnetic moment in the presence of the field 
produced by all others is given by [3, p. 444] 


^ J d' x h M • m = - J d x |h M | 

where its first variation is obtained using the reciprocity relation [3, p. 387] 

SeM = — - f d 3 x5h.M • m — - f d 3 x\iM ■ 5 m 


= — / d 3 x h m • 5m 


(2.34) 


(2.35) 


The energy expression (2.34) allow us to know what equilibrium configurations are ex- 
pected from LL update only with the demagnetizing field and to explain the second picture in 
figure 1.2. As shown, if we assume uniform magnetization along the y direction then within 
the system p* = 0. On the bottom of the system we will have negative density of charges, 
a* = (m y e y ) ■ (—e y ) < 0, while on the left and right sides a* = 0 and on the top we have 
positive density. As a consequence the demagnetizing field will point downward, opposing the 
magnetization. Since the field is nonzero, from (2.34) the energy will be positive, ejf > 0. 
We conclude this configuration is not an energy minimum. 

To be a minimum it is required to have no density of magnetic charges. Configurations as 
m circling the system tangent to its boundary give us p* = 0 and a* = 0, the demagnetizing 
field is then h m = 0 and the energy is ey = 0. However, experimentally it is found in ferro- 
magnets the formation of domains of aligned magnetic moments represented as an example 
in figure 1.2. The presence of quantum exchange interactions allows its formation. A new 
energy and magnetic field term needs to be considered, the exchange field. 
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Contribution 4: The Exchange Field 


The quantum exchange interaction [1] between closest neighborhoods magnetic moments can 
be represented by a local field derived from its energy formulation. From its quantum version 
[14, 15] it can derived the following equation: 


e E x = \ j ((Vm x ) 2 + (Vrriy ) 2 + (Vm z ) 2 ) d 3 x (2.36) 

v 

where magnetization and distances are dimensionless, and we note that the exchange constant 
A that determines the strength of this contribution is contained in the definition of l ex we 
use as a length scale (c/ p. 9). 


The first variation of the energy with the magnetization is given by 

2 


Se EX = b / _ ( Vm *) 2 d 3 x 

i = 0 

= - f 2 SJrrn ■ V5m.i d 3 x 


■' i=0 

where we expanded the squares and cancelled terms. Integrating by parts gives us 

r 2 r 2 

de E x = ~ E 5mi'V 2 mi + / SmjVmi ■ hdS 

d i=o d )=0 

= — j 4'm • V 2 md 3 .x + J ((n • V)m) • 5m dS 


(2.37) 

(2.38) 


(2.39) 


If this variation is close to equilibrium, as we shall see in next section, we will have (n • V)m 
close to zero, then we can identify the exchange magnetic field as given by 


h ex = V m 


(2.40) 


The functional form of the energy and field has a second derivative of magnetization, so the 
smother is the function m, the smaller are the derivatives and so is the field and energy. These 
smooth configurations are those which have a minimum of energy. Given that LL updates 
the magnetization in order to minimize energy, it is to expect that from a initial high energy 
state, characterized by a spatially irregular function m, that it evolve into a state where 
second derivatives are zero. 


The combination of demagnetizing field and exchange field yields the formation of magnetic 
domains. The magnetization is aligned making the second derivatives zero and p* = 0. In 
order to cancel fictitious magnetic charge density, a*, at the boundary of the domains, these 
are arranged for example as in the first picture in figure 1.2. 
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2.3 Equilibrium conditions 


From previous section, we identified the equilibrium configurations from minimums of the 
energy expressions, these configurations are those that will not evolve in time, = 0, by 
LL or LLG equations. From the variational point of view, these minimums are obtained from 
the first variation in the energies. By collecting terms we want the magnetization such that 
for any variation, <5m, the total energy will not vary, that is 

0 = Se to t = ~ (hap + hAAr + h M + V 2 m) • <5m d 3 x+ (2.41) 

Jv 

j (h • V)m • dm dS 

dV 


The variation needs to conserve the norm, that is |m + <5m| = 1. There are two ways 
to implement this statement, using the method of Lagrange Multipliers or instead simply 
introduce the statement directly into the variation [5]. By defining a small rotation about 56 
with dm = 56 x m, the equation (2.41) becomes 

0 = — j 56 ■ (m x h t 0 t)d?x + j 56 • (m x ((h • V)m )dS (2.42) 

V 8V 

where 

h tot. = h ap + h/iiv + hM + V^m (2.43) 


Since we require this condition for any <50, the terms in parenthesis in (2.42) need to be 
zero 


m x h tot = 0 (2.44) 

m x ((n • V)m) = 0 (2-45) 

These are the Brown’s conditions for equilibrium magnetization, the stability would require 
knowing the second variation of the energy [5] . The first condition implies that magnetization 
need to be aligned with the magnetic field in a equilibrium state, as expected experimentally. 
The second condition is 

m x (n • V)m = 0 (2.46) 

which can be simplified. Noting either m is zero, or (n • V)m is zero, or they are parallel. 
The first is excluded because the norm is nonzero, the last is excluded because norm cannot 
change, so the second must be true 

(h • V)m = 0 (2.47) 

This condition will be used in next Chapter when second derivatives for exchange fields need 
to be computed near the boundary of the system. 
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2.4 Landau-Lishitz-Gilbert equation and Rayleigh dissipation 


It is known from classical mechanics that friction and other dissipative forces can be approx- 
imated by the gradient of quadratic function in the velocities, the Rayleigh function [18]. 
Analogously, in micromagnetism we can introduce 


R = a o 





(2.48) 


If the coefficients a.; are positive then the function will represent a bowl. For a given rate of 
change, rh, the negative of the gradient will point in the opposite direction. 

Remembering that in m-space the magnetic field h determines the direction of evolution 
of magnetization m, we can introduce the dissipative force as a magnetic field. To conserve 
the norm, the update must be orthogonal to m, hence for the rotational term we have -mxh. 
To introduce dissipation on this term we can add a new magnetic field, h^s, that opposes h. 
Thus we define the following equation 


<9m 

~dt 


-m x (h + h dis ) 


(2.49) 


where 


h dis — R 


dm x dm v dm*\ 

ai, -m' ai -df a2 ^r) 


(2.50) 


We can simplify the expression by choosing the coefficients as ao = a\ = <22 = a thereby 
giving us the LLG equation 


dm dm 

—— = — m x h + am x — — 
dt dt 


(2.51) 


where the a coefficient is not the same one as in LL equation [11], 

From [16] we note that LL and LLG determine different magnetization dynamics, and 
only if both a are small these are equivalent. The LLG equation has a dissipation term, 
m x tjP, that affects both the rotation term and the damping dynamics. If we increase a 
both dynamics are suppressed which is what we intuitively expect. This does not happen 
in LL, since both terms are orthogonal, and a only affects the damping term. It is for this 
distinction of LLG from the LL equation that we choose to integrate the former in the next 
chapters. 
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Chapter 3 


Landau-Lifshitz- Gilbert equation 
and Finite Element Method 


Magnetization dynamics also satisfies LLG equation. In this chapter we aim to introduce 
Newton-Raphson (NR) method for LLG, a well known technique for non-linear equations, 
after we descretize time using the Implicit Mid-point Rule [19, 12, 10]. 

We review the Finite Element Method to discretize space and introduce the necessary 
equations to compute each magnetic field term. This method allows simulation of curved 
shapes better than Finite Differences. 

In [10] the discretized residuals for LLG and Poisson’s equation are derived using Backwards- 
Euler method for approximating Implicit Midpoint Rule and the space discretization is done 
with basis function focused on the vertices [21, p. 38] using vectorial notation. In this work 
instead we will use basis functions focused on each element since it allows a better bridge 
between the equation and the code that will be develop in later Chapters, we will use Implicit 
Mid-Point Rule to discretize time. Additionally we will prefer using tensor notation in our 
equations, the Chapter is self-contained. 

3.1 Statement of the problem 

The goal is to compute the time evolution of magnetization m(x) = (m x ,m y ,m z ) T from an 
initial configuration that satisfies the LLG equation with the applied field, the stray field and 
exchange field contributions but we will not use anisotropy as the first three are the ones used 
in the standard problem we will use to validate the code. The LLG equation holds at every 
x e IR 3 


dm 

(9m 

(3.1) 

dt 

= -m x h tot + a c m x — — 

at 

h tot 

= h ap + h M + V 2 m 

(3.2) 

h m 

= -v<t> 

(3.3) 


While in previous chapter the stray field was computed directly from the magnetization 
and hence only the system was considered and not its surrounding, in this chapter we will 
choose a new strategy and use the magnetic potential that satisfies Poisson’s equation for all 
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space, IR 3 . We will split it into the volume of the system, where |m| 7^ 0, defined to be fi and 
a finite surroundings, fi c . The following picture represents these subvolumes. 



Figure 3.1: Representation of d = 3 subvolume Q and a finite surroundings D c , by extending 
the latter to infinity, Q U fl c = IR 3 . The dotted lines trace the outline of the surfaces of both 
subvolumes. The point x belongs to the surface 


From the theory of micromagnetism the demagnetizing field potential is given by 


V 2 6 = 


V • m if x G hi 


0 


if x 0 


(3.4) 


Defining the two limits of a generic function /(x) as x tends toward the boundary point 
xq G <9D from inside ( int ) or outside (ext) as the values of the function fi n t,ext( x o) 


lirn /(x G fi) = /mt(x 0 ) 

x— >-xoE<9r2 

lirn /(x G fi c ) = fext(xo) 
x—>xq Ed ft 


(3.5) 


the matching conditions (MC) for the Poisson problem are 


dfpint d(f>ext 
— — = m • n 

dn dn (3.6) 

< fiint (ftext — 0 

where m • n is evaluated at boundary value xo. Additionally we require the asymptotic 
behaviour (j> — > 0 as |x| — > 00. The presence of Poissons’s equation and of a time derivative 
on the RHS of LLG is what makes this problem distinct from LL, it’s a implicit problem 
and new techniques such as the Newton-Raphson method and Finite Element Method will be 
used, but first we will introduce a time discretization known as Implicit Mid-Point Rule. 
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3.2 Implicit Mid-Point Rule and Newton-Raphson Method 


First we will write LLG showing its dependencies 


9m 

~dt 


f 



(3.7) 


where the explicit dependence on time comes from h ap . However, for our case, we assume the 
applied held is constant. Since time is a one dimensional variable, discretization of derivatives 
only requires replacement by finite differences 

9m m(f + At) - m(f) 

~m a t (3 - 8) 


The idea behind IMR is to compute / at the average magnetization and potential between 
the current configuration and the next one, its arguments have then the following substitutions 


( t + At) + t 


— tmid 


m 


4> 


m(f + At) + m(f) 


— ni mid 


4>{t + At) + 4>(t) 


4* mid 


(3.9) 

(3.10) 

(3.11) 


Observe that / depends now on the magnetization at the next instant of time. Substituting 
into Poisson’s equation as well we have for a given m(t) and associated (j){t ) the following 
equations 


m(f + At) — m(t) 
At 

^ < fimid 


f ( tmidi m-midi 4* midi 


m(t + At) 


At 



V • m mic i 


(3.12) 

(3.13) 


that implicitly determine the updated configuration m(f + At) and 4>(t + At). 

The Newton-Raphson strategy of finding the magnetization m(f + At) and potential 4>{t + 
At) that solve (3.12) and (3.13) is to build for each equation a relation such that these 
solutions are their fixed points, we name these exact solutions, m* and 4>* [17]. We define the 
current configuration as mo = m(f) and 4>q = 4>(t) and a generic configuration at t + At as 
m = m(i + At) and 0 = 4>(t + At) not necessarily satisfying these equations, whose residuals 
are 


r(m, 4 >) = 


m — mo 


At 


f ( 


m + mo cj) + 4>o m m 0 \ 


s(m, <f>) = V 


m + mo 


- V 


2 ’ 2 
2 ( 4* + fio 


At ) 


(3.14) 

(3.15) 


Observe that while m and 4> ar e variables, mo and 4> o work as parameters that we have to 
choose. For given ones, different residuals are obtained associated with the configurations, 
m and 4>, that measure how close they are to LLG and Poisson’s exact solutions. If we set 
up r = 0 and s = 0 for all x E 1R 3 , these two equations determine m* and (/>*, that are the 
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next step from mo and (f>o the LLG equation determines while simultaneously satisfying the 
Poisson’s equation. 

A practical approach to find m* and (j)* is to choose mo and (f) o and then approximate 
(3.14) and (3.15) by a first order Taylor series around a particular configuration m p and 0 p 
that we must also specify 


r(m, (j>) « r cuto// (m, (f>) = r(m p , </> p ) + . ( m _ mp ) + dv ^ a h^ ( tA^ _ 0 p ) 

s(m, (j)) « Scutoffi m, 0) = s(m p , (f> p ) + ^ • (m - m p ) + _ 0 p ) 


(3.16) 

(3.17) 


It’s important to notice the derivatives are with respect to the variables m and (j) and then 
evaluated at the point of expansion m p and <f> p . From these equations an iterative method 
can be devised. We start by solving r cuto ff(in, (j>) = 0 and Scuto//(m, 4>) = 0 with m p = mo 
and 4>p = <f>o, if the solution m and 4> are closer to m* and (f>*, we can repeat the process 
using these as the new m p and (j) p , progressively evolving m toward m* and <f> to (j>*, but once 
the residual is small enough we can stop the process, yielding two approximations, m 7 and <//, 
to the exact solutions. We then reparametrize the linear approximations with mo = m 7 and 
(f > o = (j )' , and refresh the process looking for the new step in the evolution of magnetization 
and potential. 

The solutions we will search for (3.12) and (3.13) are not m* and (j)* but approximated 
ones, ml and (j )% , since it will not be the residuals (3.16) and (3.17) we will use. Instead we 
will start by putting both in what we will later call as their weak forms and then introduce 
a spacial discretization. Its these discretized residuals that we will linearly approximate. For 
that we will require expressions for the discretized r(m, (f>) and s(m, (/>) and then compute 
their Jacobians yielding explicitly a 3 x 3 matrix for a 1 x 3 row vector for a 3 x 1 
column vector for and a scalar for In the next section we introduce the well known 
Finite Element Method and expressions for the entries for both residuals and Jacobians. 


3.3 Domain discretization and basis functions 

The spatial domain where the potential function is defined can be tesselated into tetrahedra, 
characterized by four vertices connected by edges. As an example take the cube in figure 
3.2, which has 8 vertices plus a central one. We can divide it into 12 tetrahedra of which 
one is shown, like all the others, it has the vertex 8 and the remaining three on one of the 
facets of the cube, vertices 4, 5 and 6 for this particular element. Observe these tetrahaedra 
fill completely the cube leaving no holes, and all edges connecting the vertices never cross, 
these properties characterize tetrahedral meshes for generic domains, observe also that the 12 
tetrahedra enclose the vertex 8 and the inner triangular facets are shared by two tetrahedra 
while the surface ones do not. 

Building meshes for arbitrary shapes as in figure 3.1 will require these two main properties, 
the shapes of tetrahedra might be different but ideally each element should have all six edges 
with the same length as this minimizes the error in FEM, not the case for our BCC example 
where 56>68 = 48 = 58>46 = 45. 
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Table 3.1: Connectivity matrix for figure 3.2. Each column represents a tetrahedral element and each 
row all vertices with local index a. 


(a,e) 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

0 

0 

2 

0 

0 

0 

0 

2 

2 

5 

4 

1 

3 

1 

2 

4 

4 

1 

2 

1 

6 

3 

6 

5 

3 

5 

2 

4 

6 

5 

5 

3 

3 

7 

7 

7 

6 

5 

7 

3 

8 

8 

8 

8 

8 

8 

8 

8 

8 

8 

8 

8 


In FEM [17] a mesh structure for SI U I2 C allows one to replace the problem of finding an 
exact solution for the Poisson equation, </>*(x), by the finding of an approximated one, ^(x), 
where the index h indicates that this is an approximated version of the exact solution within 
the tetrahedron element of index e. The approximation consists in considering the behaviour 
of i up to a given order, we will assume just a linear behaviour within each tetrahedrum, it 
will depend only on the potential on the four vertices, so we only have to find a finite number 
of values for 4>h, the n q vertices of the mesh. 

The Poisson’s equation which determines locally the behaviour of the exact solution is 
reformulated into its weak form to be introduced later, so that instead it determines the 
behaviour of <f>h by relating neighbouring vertices: for each vertex a linear equation is associ- 
ated, the set of these equations determine a linear system of equations whose solution is the 
potential at each vertex. 



Figure 3.2: A body-centered cube (BCC) with 9 vertices and 12 tetrahedral elements. Cube’s 
facet edges not drawn except for tetrahedron 4568. 

To build a system of equations that relate each vertex with its neighbours connected by 
edges, we need to represent the mesh by identifying each vertex and describe how they are 
connected. We define then the connectivity matrix for tetrahedra, a matrix with four rows 
and a number of columns equal to the number of tetrahedra, n e . 

First we index all vertices i G {0, 1, 2, ...,n q — 1}, in any order, just as in figure 3.2, and 
then the tetrahedron shaped elements with index e G {0, 1, ...,n e — 1}, which although not 
shown in the picture to simplify it, belong to {0, 1,2, ..., 11}. Then for each element, e, the 
associated four vertices are numbered again in any order, but by local indices, a G {0, 1, 2, 3} 
corresponding to the four rows of the matrix at the column e. The matrix entry (e, a) 
identifies uniquely with local variables one vertex of the mesh, the vertex we assigned the 
glocal index i. These are the entries in table 3.1, where we have the tetrahedron in figure 3.2 
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Table 3.2: Connectivity matrix for surface elements in figure 3.2. Each column represents triangle 
element and each row all vertices with local index a. The a indices need not coincide with table 3.1. 


(a,s) 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

1 

2 

4 

4 

1 

2 

1 

6 

3 

6 

5 

3 

5 

0 

0 

2 

0 

0 

0 

0 

2 

2 

5 

4 

1 

3 

2 

4 

6 

5 

5 

3 

3 

7 

7 

7 

6 

5 

7 


with index e = 9 since all its four vertices are in the tenth column. All other tctrahedra not 
shown in figure 3.2 are represented as all other columns. 

Another connectivity matrix for the surface vertices can be built. Keeping the same global 
indices for the vertices, we index the n s surface triangles in any order, s G {0, 1, 2, ..., n s — 1}. 
In this simple example in figure 3.2, the boundary connectivity matrix in table 3.2 has as 
many tetrahedral elements as there are surface elements, also the first three rows of table 3.1 
coincide with table 3.2, generally this is not the case and n e 3> n s . 

Once built, the connectivity matrix will allow us to convert from indices (e, a) into i or 
for example to identify easily all vertices with a = 2, as they all are in the third row. It 
works like a function, we will name it i = n(e, a) or i = n(s , a) depending on whether we 
are referring to tetrahedral elements or triangular elements. With them we can define sets 
like {(e,ot)\i = n(e,a )} which identifies all vertices, in (e, a) notation that correspond to the 
same vertex of the mesh and also gives us the elements they are in, this is important for 
the equations we derive next since calculations associated with a vertex i can be split into 
calculations on each tetrahedra that shares this vertex. 

Now that the mesh structure is well defined we can identify the four vertices of an element 
e by the four global indices in column e, we then easily pick their locations by preconstructing 
an array with three columns, for x, y and z and where each row is assigned i. With these 
four vertices locations we can obtain the linear representation of the potential or any other 
function of space within the tetrahedra. We proceed as follows for the potential case, define 
it within the tetrahedra e as 0^(x) with x G I2 e and suppose it can be written as 


<Mx) 


a e + b e x + c e y + d e z if x G 

o if x 0 n e 


(3.18) 


where a e , b e , c e and d e are unknowns to be determined that characterize the function inside the 
element. Close neighbour elements of a given vertex i. { (e, a)\i = n(e, cc)}, need to output the 
same value for </>(xj) = 4>i to give us a continuous solution. To sustain this consistency through 
out the mesh we need the values of (f>i to determine these coefficients, in local coordinates the 
four potential for a single tetrahedron are </>q, (j) f, 0| and the following system of equations 
allow us to relate these with the coefficients 


(1 x o Vo *o\ 


/ a e \ 



1 x\ y\ z\ 


b e 



1 x 2 V e 2 z 2 


c e 



V 1 y% 4/ 


\d e ) 


""oiCO 


22 


(3.19) 



where we define the matrix as X. 

By inverting the system of equations each of the coefficients is determined by the vertices 
potentials, . With the coefficients the potential inside the element, e, comes determined 
as well when we use (3.18), that is, the interior information is linearly interpolated from 
the boundary vertices’ potentials. We can invert the system of equations by introducing the 
cofactors of matrix X as 


( det(X) 

\ 


(l 

x 0 

y e o 

4) 


4oo 

CIO 

C02 

£30^ 


det(X) 


1 

x\ 

vt 

zf 


Coi 

Cll 

C21 

C31 


det(X) 


1 

x 2 

yl 

4 


C02 

C12 

C22 

C32 

V 

det(X) ) 


VI 

x 3 

yl 

4/ 


\C03 

C13 

C23 

C33/ 


(3.20) 


Then the coefficients are 


/ a e \ 


{ Coo CIO C 02 C 3 oV 


( 4 \ 

b e 

1 

Coi Cll C 21 C31 


4 

c e 

det(X) 

C 02 C 12 C 22 C 32 



W 


\co3 C 13 C23 C 33 / 


W 


Substituting these coefficients in (3.18) we get 


(3.21) 


M*) = def ^ y) [(coo + cio* + c 2 oy + c 3O z)0q+ 

(cio + Cnx + Ci 2 y + Ci 3 z)0f+ 
(C20 + C12X + C22V + C23^)</>2 + 
(C30 + Cl 3 X + c 32 y + C 33 Z)(/>|] 


Observe that 0^(x) turns out to be a linear combination of linear functions, whose coef- 
ficients are the cofactors of the matrix X. These are the basis function for the potential in 
element e, they are important from now on, so we name them 


KW = 


det(X) 


(4,0 + c a,l x + c a,2U + c a,3 z ) 


(3.23) 


and observe that det(X) = 6ST e where Q e is the signed volume [21], that is the volume which 
can be positive or negative depending on how the rows of X were organized which in turn 
depends on how we numbered the vertices locally. By having the signed volume, the order of 
numbering is irrelevant. 

With the shape functions we can approximate within element e the potential as 


3 

0(x) ~ <4( x ) = Y 4( x )^ ( 3 - 24 ) 

Q! = 0 

where its only nonzero within IT. These linear basis functions are also known as shape 
functions and have an important property, notice the factor in parenthesis in (3.23). It has 
almost the det(X) form. The a index indicate the row of X, where 04,y®,2:®) = x® was 
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replace by x. If we pick x = x® then IV® (x® ) = 1 in (3.23). But if x is any other of the three 
tetrahedron vertices locations or any linear combination of these, then the basis function is 
zero. Basis functions with these properties allow 0®(x®) = as required. We conclude 


= 5 a p 


(3.25) 


Not only we can expand scalar function on this basis, also vector functions like magneti- 
zation can be expanded as follows 


3 

m h = ^m®IV®(x) (3.26) 

O !=0 

We can now extend the argument and represent the functions for all tetrahedra by summing 
over all elements as 


n e — 1 3 

<M x ) = EE«W (3.27) 

e=0 a=0 
n e — 1 3 

m /t( x ) = Y i^ m a^a( x ) ( 3 - 28 ) 

e=0 CK— 0 


However care must be taken when we evaluate at a given x, (it will be useful in section 4.3) 
since all basis function incident on i, {IV®(x)|i = n(e, a)} yield the value of 1. Defining the 
number of these functions as Ni nc = #{(e,a)\ i = n(e,a)}, the potential will be = 

Ni nc x (j>i. To prevent this N mc factor when we use either (3.27) or (3.28), we will assume that 
the domain for each basis function incident on i does not superimpose and only one of the 
basis functions is actually 1 on Xj. Therefore we have 0^(xj) = (f){ and m^(xj) = m*. Notice 
this assumption is irrelevant when computing integrals on IR 3 that involve (3.27) or (3.28). 

Finally the basis has the following important property as well, for a given tetrahedral 
element the basis is not orthogonal, instead satisfies the following relation [21] 


J {N^ k (Nf) l (NI) m (NI) n d 3 x 

n e 


k\l\m\n\ 

(k + 1 + m + n + 3)! 


6H e 


(3.29) 


for surface triangular elements N t ® has the same values that a triangular basis function would, 
the latter satisfies 

/ i-l/ljT)! 

(N$f(N;)‘(NirdS = (fe + ; + m+2) , 2A» (3.30) 

Both will be useful when evaluating residual integrals in next sections. 


3.4 Weak form of Poisson’s residual 

Since within the element e, c j>h is a linear function, when we substitute into (3.4) its Laplacian 
is zero. Poisson’s equation is a strong statement since it excludes solutions of the form (3.27), 
but multiplying by a basis function, integrating we get 
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V • mN^(x)d 3 x 


(3.31) 


J V 2 </>l\T(x)d 3 x = J 

n e n e 

and while before we had the Poisson’s equation (3.4) relating the local curvature of the 
function cj) with the divergence of the magnetization and this determine the solution, with 
this integral equation (3.31) we have on its LHS the sum of second derivatives weighted by 
the shape function, equal to the RHS where we have the sum of divergences of magnetization 
also weighted by the shape function on element H e . This statement is less restrictive because 
takes into account information from the entire element, it determines a larger set of <f> solution 
for (3.31) [10], 

However to allow a linear solution we still have to go further. Integrating by parts, we 
can transfer one of the gradients to AT(x) at the price of a minus sign and a boundary term, 
rearranging we define the Poisson’s residual in its weak formulation 


= J V • mJV'(x)d 3 i + Jv<p- VAT(x)d 3 x - J A £(x)V0 • MS (3.32) 


When we set it equal to zero we obtain the weak form of Poisson’s equation. However our 
goal is to implement Newton-Raphson, hence its the residual that we will work with. We can 
transform it into a set of linear expressions by substituting fa and m^, this is possible since 
the second integral now isn’t zero, instead of the Laplacian we have the gradient of the basis 
function which is given by 


VAT(x) = 


2 ^ 1 , 0 ! 
{x) = Mi(Xj C2 ’“ 

\C3,a 


(3.33) 


From (3.32), the residual at vertex i can now be computed, summing contributions from 
elements that share this vertex we have 


^ = E ( 3 - 34 ) 

{(e,a)| i=n(e,a)} 

Notice that before discretization we had for every function ^>(x) we substitute in (3.15) a 
residual value for every x £ IR 3 , now after discretization we have it for each vertex of the 
mesh, Si . Arranging it into a vector we get 



which is the discreet version of the continuous function s(x) in (3.32), information about 
its behaviour is lost but now we only need to handle an n q entry vector, which can better 
approximate it’s continuous counterpart if we reduce the element size and thereby increase 
the number of n q of entries in (3.35). 


25 



Each one of the three integrals in (3.32) is a contribution for the i = n(e, a ) entry of a 
column vector, they behave differently depending of the the location of vertex, x,;. We will 
analyse them starting with the source term (src), followed by the stiffness term ( stf ) and 
finally the surface integral ( d ) where Neumann boundary conditions and matching conditions 
enter. The sum total is given by 



(3.36) 


However, if we suppose we have instead a Poisson’s equation with Dirishlet boundary con- 
ditions on the boundary of the system, dQ, then there would not be a surface integral and 
residual contribution i to take into account [21, p. 92, 108]. 


1. The Source term 


We start by approximating the divergence using the discretized magnetization 

3 

V • m ps V • m h = • ViV* (3.37) 

a=0 

the source term integral in (3.36) is given by 

3 

V • m/jIV|(x)<i 3 x = — ^ • ViVj (3.38) 

4 0=0 



where we used the fact f Qe N^(x)d 3 x = from (3.29). This contribution only depend on 
e. To sum all contributions for the vertex i is to sum factors for each e that have i on its 
vertices. 


src 


E 

{(e,a)| i=n(e,a)} 


S <t>,a ~ 


E 


y fie 


{e|3a(«=n(e,a))} /3=0 


m 


VJV| 


if i G fi U fi c 


(3.39) 


Observe that for i e all terms of the double sum are nonzero since this is completely 

within the magnetized material. If i £ dQ then terms associated with elements outside the 
magnetized region will be zero. For i within the vacuum region, then the source residual term 
is zero, we have a Laplace’s equation. 

Substituting now IMR, we have 

3 fi e 

with IMR: 47 = y ( m /3 + m ) 3 ,o) -VJV| (3.40) 

{e|3a(2=n(e,a:))} fi=0 

There are n q residuals like these, as many as the number of vertices in the magnetized region 
plus the surrounding vacuum, with them we can build a vector, s^ rc , which will be summed 
to the other two integral terms we analyse next. 
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2. The Stiffness Term 


Introducing the descretized potential from (3.24) into the second integral of (3.32) we get 

/■ 3 

45 = / V0 h ■ ViV^(x)d 3 .T = L! e £ (ffp\7Np • VIV* (3.41) 

/ 3=0 

Summing now all contribution from all elements that have the vertex i we have 

<i"= E = E X>JVJVJ • VAr»ST if i£f!Uf! c 

{(e,a)|i=n(e,a)} {(e,a)|i=n(e,a;)} /3=0 

(3.42) 

1 3 

with IMR: s£/ f = - £ J](^ + ^ >0 )ViV^ ViV^ e (3.43) 

{(e,a)|i=n(e,a)} /3=0 

Many terms of this double summation in (3.42) will have <j)p corresponding to the same 
c t>j , where j is a global index of a vertex connected to i. Lumping all its coefficients into a 
single constant, we define it as Kij. 

K i:j = VN e p-X7N^Q e (3.44) 

{(e,a,/3)|i=n(e,a)A.j=n(e,/3)} 

We see then, s s £ l f as a sum of products of each K l j with the corresponding 0p By going 
further we extend this sum by adding the multiplications of K tJ = 0 with 0j for all other j 
not connected to i. This shows that we can express s s ^ as a dot product of a row vector 
with Ki j as entries with a column vector 0. Since all other indices i could have been chosen, 
we have for each one a dot product like this, all these calculations can be arranged in a column 
vector s% tff which is then given by the product of the matrix K with the vector 0 

»f" = ^(* + *0 ( 3 - 45 ) 

For a given row i of K the only nonzero entries are those whose columns that are associated to 
vertices connected to i including the diagonal entries that correspond to ^{( e a)\i=n(e «)} VJV-. 
V1V®. Since most of the entries of K will be zero, we designate K as a sparse matrix. 

The n 9 xl vector 0 is multiplied by the n q x n q sparse matrix K to yield the second piece 
of s for (3.35). 

3. The Surface term 

The surface integration is made along the four triangle facets of a tetrahedron, (9fl e 

45= / Ar a(x)V0 fe • MS (3.46) 

an e 
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The surface contribution for a vertex i is obtained like the two previous integral terms by 
adding the contributions of all tetrahedra that have i as one of its vertices, we write 

4 = E / KWVt-ndS (3.47) 

r / \r • ( w J d£l e 

{(e,a)| i=n(e,a)\ 

The result of such a calculation depends on where the vertex i is located. We have four 
cases to consider. 

If i E £4 \ <912 or i G fl c the surface is composed of triangles that form a ’’shell- like” 
surface that is integrated once with n pointing outwards, analogous to i = 8 in figure 3.2. 
Remembering the basis functions are equal to 1 at the core of this shell, x$, all basis function 
will be zero on its surface. Integrations along triangular facets inside the shell will be made 
twice with h pointing in opposite directions, since X7(f> ■ n is continuous within or outside f2, 
these terms cancel in pairs. Hence s ® ■ will be zero within these two subvolumes. However if 
i e <912, some of these facets belong to (912, as an example consider the following figure where 
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Figure 3.3: Four elements with shared facets belonging to the tessellated surface <9f2 (shaded). 


we represent only four of the elements that will enclose the central vertex 2 = 1 that belongs 
to the surface <9f2. Two terms of the summation in (3.47) for this case are 

4=i = [ • niVr°(x)dS + f V <p ext • (-n)ATr 1 (x)d5 + 

7Ai 23 2Al23 

= f m • nN^ =0 (x)dS + ... 

J A 123 

where n(e = 0 , a) = n(e = l,/3) = 1 hence both basis function for different elements will 
give the same output if evaluated on the same triangle, A 123 , by (3.6) the derivatives on <912 
are discontinuous, therefore the integrals do not cancel. Generically we will have to sum all 
integrations along each surface triangle in <912 that has i as one of its vertices, we then write 
the nonzero part of (3.46) as 

4 = J As m h ■ nN^dS (3.50) 


(3.48) 

(3.49) 
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Summing for all contribution for vertex i 




mft • hN^(x)dS 


2 

Y m /3 ' x 

0=0 


12 


A s 



if a/(3 
if a = /3 


if n(s, a) 6 90 


(3.51) 


with IMR: 


2 

^(m^ + m^ 0 ) -h s x 
/3=o 



if a / /3 
if a = P 


(3.52) 


where observe we substituted = Yp= o m^IVg, that is all (e, (5) vertices also belong to 90 
and we used the integral (3.30). 

The final case occurs when vertex i belongs to the vacuum domain surface 90 c , figure 3.1 
the vertices connected to i will be either within 0 C or on the surface, 90 c . This case is similar 
to figure 3.3 if we suppose there are no tetrahedra in 0 C and we substitute O — > O c , 90 — > 90 c 
in the figure. Here we will only have surface integrals on the internal side, to compute them 
we need to specify the derivatives, V0 • h. Since we require our solution to decay as |x| — > oo, 
if we choose 90 c very far from O we could set the values of the derivative to zero, but instead 
we can truncate this vacuum region closer to 0 and introduce the expected behaviour of the 
potentials derivative, for d = 3 we have [21, p. 184] 


(V</> • h)(x) 


1 

x 


</>(x) for x e 90 c 


(3.53) 


Assuming |x| is constant along the triangular element of tetrahedron e that also belongs 
to 9f l c we suppose its the mid-point of triangle and call it x^ id 


q d = 


£ 


{ (e,a) \i=n(e,a)} 


[ ^-fa( X )N^dS 

ldfl e l x l 


£ 1 

{(e,a)|i=n(e,a)} 1 mtdl 0=0 J dU 


N e a N e p dS 


£ 


{ (e,a) \i=n(e,a)} 


X 


m.id\ ^_ 0 I |A e if a = /3 


if i e dn c (3.54) 


with IMR: = 


£ 


{(e,a)| i=n(e,a)} 


X' 


e I 

'midi «_ q 


s i J.S 


+ 


^A e if a / /3 


3,0 1 


24 

^A e if a = f3 


(3.55) 


where we substituted the discretized potential fa, which collapse into an expansion on just 
three basis functions instead of four since we are computing </>/,. on the surface of a tetrahedra. 
The areas A e are for the triangular elements belonging to 9fl c and the element e. From the 
n q terms j we build the final part of (3.35). 
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3.4.1 Jacobian of Poisson residual 


Now that we have expressions for each part of Poisson’s residual the derivatives can now be 
computed. Taking into account the signs in (3.36) and the expression that only depend on 
m, (3.38) with IMR and (3.52), we have the Jacobian 


with IMR: Q e ai 


ds e a 

9m( 


Q e 

—VAT - n s x 



if a 7 
if a = 7 


(3.56) 


as we observer in (3.55), the second term only appears if (e,a) and (e, 7 ) both belong to d£l 
This first Jacobian is a 3 x 1 vector, but as we shall see in Chap 5 it is more suitable having 
it as a row vector, hence we transpose (3.56), it is associated with the indices i = n(e,a ) 
and j = ra(e, 7 ) that identify its location within a larger Jacobian matrix 5 n q x 5 n q to be 
introduced later. This term only depends on j hence its the same for all i in the matrix 
column j. 

The second Poisson’s Jacobian is obtained applying to s^ a , one of the parts (3.36), 
noticing s® is linear function with respect to , its coefficients are the derivatives, therefore 
from the stiffness term (3.43) and if BC (3.53) are used in (3.55) we have 


if a 7 ^ 7 and n(a, 7 ) £ dQ c 
if a = 7 and n(a, 7 ) £ dQ c 
if n(a, 7 ) 0 <9fl c 

(3.57) 

where the derivative turned (j) s p 0 and <^j| 0 into zero, so there is no contribution. This term is 
a contribution to the entry of K matrix obtained previously, but now it includes as a second 
term the asymptotic boundary conditions. 


-K e 
2 " 7 


+ O e = 


94 


d(f>i 


Ivk-vn^ + I- 


1 


X < 


mid I 


J-A e 

24^ 

J_A e 

12^ 


3.5 Weak form for LLG equation 


In an analogous way to what have been done with Poisson’s residual, the weak form of LLG 
residual is obtained following the same procedure. We multiply LLG equation (3.1) with a 
basis function, integrate and rearranging we define 





Am h 

At 


+ m/, x 



V <4 + a c 


Am fe 

At 



i4(x)d 3 .x 


(3.58) 


It is more complicated than the Poisson’s residual, but we will separated it into three parts, 
the time derivative, the applied field together with the stray field and damping, and finally 
the exchange term. Some of the integrals involved will be approximated using the quadrature 
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from [10]. For a generic function /(x), its integral in an element of volume fl e is approximated 


[ f(*)d 3 X « ~F/( X a) ( 3 - 59 ) 

Once we have r® and its Jacobians we can assemble the residuals at i as the following sum 

r i = Y r « ( 3 - 60 ) 

{(e,a)\i=n(e,a)} 

and an entry of the 3 x 3 or 3 x 1 Jacobians as 


drj 
dm j 

dr, 

d<l>j 


E 


dr e 

a 


— dm.% 
{(e,a,/3)|i=n(e,a)Aj=n(e,/J)} P 

dr e 

a 

— dfin 

{(e,a,/3)\i=n(e,a)Aj=n(e,/3)} P 


E 


(3.61) 

(3.62) 


The computation of such sums will be one of the subjects for Chapter 5, for now we will 
just derive their factors from the expression (3.58). 


1. The time derivative 

The time derivative residual of (3.58) after discretization by finite differences is given by 


r e,A _ 


m h - m hfi Me , 3 


At 


N e d 6 x 


Q e ^ m /i ( x ^) - m fe; o(x^) 
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At 


KiO 


j=o 

n e m h {-x. e ) - nu o(x® 


4 At 

where the integral was approximated using (3.59). The Jacobian is simply given by 


(3.63) 


dr' Q 


e,A 


A e = 

Q ’ a dm% 


= I 


3x3 


i n e 

At 4 


(3.64) 


using the fact = ^ 3 x 3 and (3.29). This term fills the diagonals of the 3x3 blocks 

associated with ( i , i) with a multiple of the constant term -^j yp . That multiple is the number 
of connections of a given vertex with its neighbours. 


2. Applied field, Stray field and Damping Jacobian, 

The residual from LLG is 

hap + h a p ;0 _ / <\>h + <%,Q 

2 I 2 


mV, — m 


— CL c 


h, 0 


At 


N e a <Px 


(3.65) 


r e ’ D = 


m, 


+ m io 
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Using the integral approximation (3.59) we get 


r e ’ D = 


m a + m a,0 


h e + h e 

±a -ap,a 1 i 


ap,a, 0 1 

2 “ 

0=0 


+ 0^ o )VlV|-a, 


m a - m a,0 


At 


(3.66) 

Observe that the magnetization is evaluated at x® while the gradient of the potential remains 
unaffected by the integration since it does not depend on x. This residual term has a depen- 
dence on magnetization uih and potential field 4>h, two derivatives have to be computed. The 
first one is given by 


_ dr% _ d f + 

a ’ 7 9m( J n e V 2 ) 

dv e 

and we will compute by components. 


hap + h 


ap,0 


- V 



— a 


m h - m h , o 


At 


md 3 x 


(3.67) 


Using Einstein summation for lower indices, and knowing [a x b]j = eij k o-jb k , we observe 
that the rate of change of the i £ {x, y, zj component of r® at vertex n(e, a ) with respect to 
the l £ {x, y, z } component of magnetization at vertex n(e, 7 ) is given by 


d r Qti _ d f m h j + m h ,jf,0 f hap,k hap,k, 0 d (f>h + 4 >h, o rnh,k ~ rn h ,k,o \ r ,3 

dw n .i drru h i J n e 6ljk 2 \ 2 dx k 2 ° At ) ° X 

(3.68) 

where we dropped the index e. From expansions for m/j and <f>h, the derivatives are 


dm h ,k 

drn- h i 

94>h _ 
dx k 


= hiN* 



(3.69) 


using the product rule on (3.68) and substituting (3.69) 


d^a,i f N a AT / h apk + h a p } kfi 0/3 T 0/3, 0 ^h,k m , h,k,o\ 

drruj €ink 2 \ 2 2 dx k ° At ) 

7 371 /jj + Tflhjfi ^3 

~ a€ijl ~At 2 ‘ T 


(3.70) 

(3.71) 


Now switching indices €ij n = —e in] and renaming the index j as k we can factor out eu k in 
both terms, combining them we have 


dr a p 

dm n ^ n 


einkKK ' - 


hap,k T h a p^ k fl 


+ 


+ ‘ 


to 9N% 


dx k 


— a 


frih,k, 0 

At 


<Tx 


(3.72) 


In our case the integrand in (3.72) will be evaluated at each x® of tetrahedron e giving 


dr e ■ 

a.i 


dm: 




-cuk Y N a(*p) N j( x i3) v k(xp) 
/ 3=0 


n e 

T 


(3.73) 


where v k with k £ {x,y,z} corresponds to the three terms in parenthesis in (3.72). The 
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summation terms are only nonzero when a = (3 = 7 therefore we have 


dm a,l 


~ —tilk 


hap,k hap,k , 0 


+ 0/3,0 _ rn ht k , 0 

4 ac A t 


(3.74) 


This calculation can be represented in vectorial form by observing first what the operation 
—6uk does to the factor in parenthesis in (3.74), that we define as Vk 
Construction the following 3x3 matrix block 


0 -v z Vy 

v z 0 —Vo 


(3.75) 


The specification of an entry in this skew matrix block is made with i,l £ {x,y,z} in (3.74). 
The tensor operation — euk transforms 17 into an entry of this matrix block 3x3, and the 
block in turn is associated with (n(e,a),n(e,a)). As we shall see, these blocks are diagonal 
blocks in the larger Jacobian matrix, meaning the residual at vertex n(e, a) only varies when 
the magnetization at this vertex varies as well and does not depend on neighbour vertices. 

In vectorial notation, the (3.75) block is given by 


ne _ dr a _ & „ KpiKt) + h ap,0« 

Da > a ~ d^ a ~ T skew 4 


W^/3 + ^/3,0 V¥ e «c / e 

+ 2_j 7 ViV /3 “ X7 m ^,o(x 0 


(3.76) 

where the skew operator maps the three entries of v into the 3x3 skew-symmetric matrix 
given in (3.75). 

Using similar techniques to the previous derivative we have 


pe ex _ 

a ' p ~ 0o<:. ~ doi 


x V («±«£) jv«(x»)A = 


-j p ^ 

= - A£(x e )(m| + m^ 0 ) x 2 —fv N^x 

4 /3=0 d( Pi 

= -7 f K(x e )( m h + m h,o) X 

1 3 ne 

~ - 4 + m k,oK)) X VA 7 

p=0 

= -^ (m a + m a,° )xV ^ 7 e 


(3.77) 


' p will contribute for a 3 x 1 column vector dependent on both i = n(e, a) and j = n(e, 7). 


3. Exchange Jacobian, 

To compute the Jacobian of the exchange term in equation (3.58) we will treat one of the 
residual components, i. Using summation convention on i,j,k and p, q, all belonging to 


33 



{x, y, zj and dropping the upper index e, we start integrating by parts and substitute rip, 


r e • = 
' a, i 


I it 


m x v m) 


. N n d 3 x = — 


+ 


f d(m h j N a ) dm h k 3 
/ £ijk ^ q d x 

VXp OXp 

f AT e dm h , k 

I o Tiqdo 

!dn e ° X q 


(3.78) 


If we assume Brown’s second equilibrium condition (2.45), the surface integral integrand is 
zero 


dm 




h,k 


dm 


dXn 


— tijk'M'hJ ' 


h,k 


dh 


m/, x 


dm-h 

dh 


= 0 


(3.79) 


Using the derivative of a product rule on the first integral of (3.78) we have 


dN a dm h , k 3 f 

tijk'W'hJ o o d X I 6ij]^N c 

OXn OXr) / oe 


dm h j dm h , k 3 


dx p dx p 


d 6 x 


(3.80) 


The second integral is of the form a x a = 0. One of the components k £ {x, y , z} of the 
expansion rri/j is given by 


mh,k = ^ fny,kNp = mp jk Np 
(3=o 

After substituting (3.81) into the first integral in (3.80) 

f dN a dN VAT ^ 3 

r a ,i = - eijkmppmp^— — Npd x 

Jfl e OXp OXp 


(3.81) 


(3.82) 


Except for Np no other factor depend on x hence using either (3.29) or the integral using 
(3.59) we have after recovering the sums 


Trvi — 


I E m /3 x m v I • V1V,,— = “ X m v )iVN a • VN v — (3.83) 


>/=0 \p=o 


i}= 0 


where we make the important observation that entries from the Poisson’s stiffness matrix 
(3.44) are also used in the construction of the exchange residual and for a given element 
is 


m !ot = E m /3 

(3=0 


(3.84) 


Hence renaming the dummy variable 77 as (3 and joining the three i components into a vector 

(3.85) 


r " = ~4 


4 m toi x ^E m }K,/3 j 


with IMR: r e a = -- A 
4 


1 mf f + mf 


Hot w ‘■‘■Hot , 0 


3 m| + m| j(J 


E 

.(3=0 


I\, 


a, (3 


(3.86) 


where we included the in K a . 
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To compute the Jacobian we take derivative — of (3.82), using the product rule 


dr a:i 
dm. 


— ^ij k / N Q 


7 ,n 


Substituting = <5/3 7 <5,- n and = S rn 5 kn , replacing e un = and renaming j as A: 

dr a ,i 


3 derivative 

x o' 

d N a 

drnpj 

dx p dx p 

<9ra 7?n 


dm. 


rj,k 


dm. 


7,n 




(3.87) 


r Ar aiv Q a/v„ . . , l3 

o — ^ink I \pT"q,kO p'y '^ r l(3,k^p'y\ d x 

cra 7?n Jn e c,x p c,x p 

Changing the basis indices applying the Kronecker delta’s we have after rearranging 


(3.88) 


dr at i 
dm.. 


— ^ink 


7,n 


f dNg 

In* dx p L 


A7 dN v A7 <91V 7 

'Wl'rpkN'y „ TTl/3,kNf3 


— £ ink 


dx p 
dN a dN v 


(lX 


m " n,kN ' r ~dx~~dx m P’ kN l 3 


9x r . 


dN a dN 7 
dxp dx p 


d s x 


(3.89) 

(3.90) 


On the both integrals only N ’® and No are dependent on x, its integration gives 4 ' factor. 
We conclude that the component i G {x, y, zj of the 3x3 Jacobian block associated with the 
indices (n(e, a), n(e, 7 )) is given by 


dr. 


a,t 


dm. 


^-ink 


7,n 


^ 77 , k 


dN a dNrj 
dx p dx p 


tot,k 


dN a dNo! 
dxp dx.p j 


n e 

T 


(3.91) 


We can rewrite this equation using the stiffness matrix as 


dr ® 


f)m e 

U /, 6 7,72 


^ink 

'~ 4 ~ 


^ ^ 0,k^ a,0 ! ^ l.ot .k ^o.~ 

0=0 


(3.92) 


Each component of the vector in parenthesis is mapped into the skew matrix by — €i nk , we 
write then the 3x3 block associated with the indices i = n(e, a ) and j = n(e, f3) as 


E F ar . = skew 


with IMR: Et = skew 

Lt , 1 


E> 

P = 0 


K 


a,0 


— Ill 


K e 

e ± '-ag 

tot A 


m 0 + m A0 K a,0 _ m tot + m tol,0 K e an 

2 4 2 4 

0=0 


(3.93) 

(3.94) 


where the Poisson’s stiffness entries are again used in this term. 


3.6 Conclusion 

In this Chaper we laid out the plan to solve LLG and Poisson’s equation. We first discretized 
time with increment At, using Implicit Mid-Point Rule: derivatives with respect to time were 
replaced by finite differences and function dependencies on time, magnetization or poten- 
tial were replaced by the average between two successive instants of time. Then space was 
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discretized. Permeating space with a d = 3 mesh with n q vertices that defined n e tetrahe- 
dral elements, allowed us to introduce four basis functions, IV®, for each one, and for surface 
triangular elements s, the basis functions N£. The potential and magnetization were then 
approximated by expanding them in this nonorthogonal basis as in (3.27) and (3.28), which 
upon substitution the strong equations LLG and Poisson, showed us these cannot be solutions. 

In response, we multiplied them by a basis function and integrated by parts on Ci e , we 
then obtained their form without second derivatives and only gradients, allowing now linear 
solutions. 

In order to solve these weak equations for each vertex i we introduced Newton-Raphson 
method which for a given initial configuration linear approximates the residuals of both equa- 
tions for each i yielding a linear system of equations where each equation is associated with 
a vertex of the mesh. Iterative search for zero residuals with this linear system progressively 
gives us better solutions for the configuration for t + At that satisfy both weak forms. 

This system of equations is formed by the Jacobians of residuals for both equations, we 
derived explicit expressions for Poisson’s residual and for LLG residual, the time derivative 
Jacobian, the applied field, stray field and damping, and finally the exchange Jacobians. 
These depend on whether vertex i is in Q or or the boundaries. The construction of these 
Jacobians depend on having first a mesh for fl U Q c . We observe that Poisson’s equation was 
the only term requiring a mesh in I7 C since we have asymptotic behavior to be satisfied. 

What if we could avoid meshing fl c ? Then we would solve the LLG equation and Poisson’s 
equation focusing only on the system. In the next chapter this is what we will do introduc- 
ing how equations from the Boundary Element Method help to solve the FEM problem for 
Poisson’s equation. 
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Chapter 4 


Finite Element Method and 
Boundary Element Method for 
Poisson’s Equation 

4.1 Introduction 


This chapter is about the Poisson’s problem. The goal is to reformulate it in such a way 
as to replace the asymptotic boundary conditions by Dirichlet boundary conditions on the 
boundary of the system, dfl. As a consequence we will not need the system surrounding, 0 C , 
in computing the evolution of the magnetization and potential. 

The central idea is to introduce a new potential, u, that can be computed with FEM 
only on Q, and whose boundary values can be mapped into the (f> at the boundary as well, a 
strategy used in the Boundary Element Method, [29, 30]. By knowing (j> on dQ, we only have 
to solve Poisson’s equation in Q. 


4.2 A new potential, u 


We start by splitting as (j> = u + v and substitute in both (3.4) and the matching conditions 
(3.6) we get two new equations. Defining conditions on the potential u [30], conditions on v 
are obtained. Suppose that 


V 2 u = V • m if x S O 
u = 0 if x S O c 


with Neumann boundary condition 


duint 

dn 


= m ■ n 


if x € dil 


(4.1) 


(4.2) 
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As a consequence of (3.6), the asymptotic behaviour for </> and the assumption on u above, a 
few lines of algebra will show that the potential v has to satisfy 

V 2 u = 0 for xGllUff (4-3) 

d v int dv ex t . 

— — = 0 for x G <912 (4.4) 

on on 

Vint - Vext = ~ U int for X G 912 (4.5) 

v(x) — > 0 as |x| — > oo (4-6) 

The goal is to establish a relation between </>mt( x ) and rq n t(x) on the boundary <912. The 
strategy will be to build an integral equation to obtain Uj nt (x) from % f (x) using (4.5) and 
then substitute into 4>i n t = Ui n t + Vi n t- 


The first step is to relate all values of v(x) for x G 12 with all v mt (x) using an integral 
equation and for that we will use the ’’free-space Green function”. This Green function will 
satisfy the fundamental equation [3, p. 251] 

V ,2 G(x, x') = -<5(x - x') for x, x' G 1R 3 (4.7) 


and asymptotic conditions G(x, x 7 ) — > 0 as |x — x'\ — > oo. The solution of this PDE can be 
obtained using the divergence theorem, for d = 3 we have [3, p. 252] 

G(x, x) = T — for x,x' G IR 3 (4.8) 

Att |x — x / | 

We now define the integral equation for (4.3) similarly as done in FEM, (3.32). While in 
FEM, the weight functions were the shape functions, JV®(x), in BEM, the weight function is 
the free-space Green function above. From (4.3) we multiply by (4.8) and integrate 



V 2 v{x')G{x,x')d 3 x' = 0 


In an analogous way to FEM we now use integration by parts (3.32) giving us 


(4.9) 


f V / u(x / ) • V ; G(x, x!)d 3 x' + f G(x, x / )V / u(x / ) • h'dS' = 0 (4.10) 

Jn Jan 


But now in BEM, we proceed by using integration by parts again to introduce a new surface 
integral 


n(x / )V /2 G(x, x)d 3 x' — f v(x')'V'G(x,x) ■ ndS' + [ G(x,x)'V r v(x) ■ h'dS' = 0 (4.11) 

: Jan Jan 


It is the first integral in (4.11) that will require more attention. If x G 11 \ <912 then using (4.7) 
and the Dirac Delta function property 



/(x) if x G 12 \ <912 


(4.12) 
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the first integral yields v (x) and the expression 

v(x) = — f v(x')'V'G(x, x ; ) • n'dS' + f G(x,x')'V'v(x') -n'dS' if x 6 fl\tlll or 12 c (4.13) 
Jan Jan 

This expression uses information of v(x) and Vu(x) on the boundary for ( int ) or (ext) values 
to compute v(x) within the system 12 or outside. It is possible from surface information to 
evaluate interior or exterior points because these two surface integrals are in fact a reformu- 
lation of the volume integral in (4.11), and this integral is of the form (4.12). However what 
we want is to relate the boundary values Vi n t(x) to substitute into MC, (4.5). Since on the 
boundary the Delta function property cannot be used, this suggests that we extend the 12 
region to include x, let it be a semi-sphere centered at x and radius e such as illustrated in 
figure 4.1-1. Following [29, p. 49-51], we first relate the Vi n t(x) values by evaluating (4.11) at 



Figure 4.1: Locally at x of figure 3.1 from previous Chapter. Top figure (I) represents the 
extension of 12 to include the boundary point x; The lower picture (II) represents the extension 
of 12 c to include it. In both we show the case where the surface <912 is smooth, so locally the 
surface looks like a plane and the semi-spheres with radius e are half-spheres, both determine 
two complementary solid angles of 27 r. 


x G c212. Then we extend 12 to include x using a semi-sphere. Now (4.12) can be used on the 
first integral. We define the semi-sphere surface as <912 e and the remaining surface with a hole 
as <912_ e , see figure 4.1-1. Using the Delta function property on the first integral in (4.11) we 
get 

%t(x) = - J Vint(x)V'G(x,x) ■ h'dS' + J G(x,x')Vv int (x') ■ n'dS' (4.14) 

9f! e U(9r2_ e 9f2 e u9f2_ e 

Both integrals can be split into two integrals: on c212 e and c212_ e . We will compute first the 
one on the semi-spherical region, <912 e . Defining R = |x — x 7 | and using spherical coordinates 
we have e# parallel to the outward normal to the surface of 12 as in figure (4.1), and from 
V' = V x /_ x = e R ^ + e e + <§0 Rs l n ( 9) ^ we have n- V x /_ x = ^ . The first integral is then 
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given by 


lan e 


n int (x')V / G(x,x / ) • ndS' 




1 

AnR 2 

#0 

4ttR 2 



Vint{*) 


(4.15) 


H(x) is the surface area of the semi-sphere that extended Q into D c . Observe that is the 
solid angle, figure 4.1 shows an example for a smooth surface, the solid angle is then 2n. Since 
(4.15) does not depend on R then it is valid for any R, therefore 


lim f n int (x')V / G(x,x / ) • n dS' = - 7 ' n (x)u int (x) 

Jan. 


(4.16) 


where y n (x) = 

For the second integral along dtt e we have 


[ G(x,x')V' Vi n t(x.') ■ n' dS' = / 

Jan t Jan e 


1 1 
47T lx 7 — X 


Vv int (x') ■ h'dS' 


= W If Vv int {x!)-h'dS' 
4 7T R J d n 


1 an, 

R , . „ 

— \7v in t{]t) ■ n 


(4.17) 


Taking the limit we have 


lim [ G(x, x.^'V'vint (x 7 ) • h'dS' = 0 (4-18) 

Jan e 

Since both the third and forth integrals are well behaved [29, p. 51] then we have the 
integration along — > dQ. We now reformulate (4.14) with (4.16) and (4.18) as 

%t(x) (1 — 7 n(x)) = — [ v int (x')VG(x, x') • n dS' + f G(x, x / )V / n int (x / ) • n' dS' (4.19) 

Jan Jan 

By comparing the case where x is inside D in (4.13) with this expression, we can qualitatively 
interpret that when the Delta function is on the boundary ’’part of it” is outside fi, so ’’only 
part” of Vi n t(x) is picked when computing the two surface integrals in (4.19) while in (4.13) 
all of ni n t(x) is the output. How much the Delta function is left out of D is measured by the 
area of the semi-sphere outside the original D with respect to the total area of the unit sphere, 
this is what 7 n( x ) measures in (4.19). 

Now we need to obtain an analogous relation for v ex t(x), we will start with D c still finite. 
Noticing that it is not a Jordan region [31, p. 996] the divergence theorem cannot be applied, 
so we divide it into two Jordan regions: Df and Dg, as in figure 4.2, we can now apply the 
divergence theorem and also integration by parts. Now choosing the same point x on the 
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Figure 4.2: Two Jordan regions: O'j, The vectors are normal to the surface and the 
sequence 1,2,3 shows part of the outline of boundary 


surface now defined as diV\ we extend to contain it as in figure 4.1-II, the first integral 
in (4.11) is well defined. We then obtain for fl); an analogous expression of (4.14) where the 
boundary is again split into the semi-sphere e and the rest as _ e . 


v ex tW = 


1 


V ext {x')\7'G{x, X 7 ) • fl' ext dS' + 


G(x,x / )V't> e xt(x') • fi' ext dS' 


dQ\ _ e 


dQf e UdQf _ e 


(4.20) 

The integral relation for is also °f the form (4.11) but since x is on the first integral is 
zero, rearranging we have 


0 = - 


v ext (x')V'G(x, x) ■ h' ext dS' + 


G(x, x')V' v ex t{x!) • n' ext dS' (4.21) 


Summing (4.20) with (4.21), integrations along the shared facet cancels as n' ext points in 
opposite directions, figure 4.2, only the integrals along the internal facet, designated as <90i U 
<90 2 = 90, and the outer facet of <90 c are left. 

Now extending the outer boundary of <90 c up to infinity, since v(x) — > 0 and Vr(x) 0, 
then only integrals along 90 = 9 O e U 90_ e remain 


lirn v ext (x) = - f v ex t(x)'V'G(x,x)-h' ext dS , + f G(x,x')Vv ext (x')-h' ext dS' 

dw^-oo J d n £ jon_ e Jdfi £ uan_ e 

(4.22) 

Once more we have four integrals similar to those already computed for 0, (4.16) and (4.18), 
the first will give us 

lim [ v ex t(x')VG(x,x') ■ h' ext dS' = -^f u v ex t{x) (4.23) 

Jan. 

where 47T7 u is the solid angle complementary to 4vr7n obtained previously, both summed they 
complete the surface area of unit sphere centered at x, implying 7u + 7n = 1- 
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The second integral along <9I2 e in (4.21) is analogous to (4.18). Finally both integrals 
along <9I7_ e reduce to integrals along df7 as e — > 0. Observing that h' ext = — h for the region 
of dQ.\ _ e coincident with the original <917 as shown by comparing figures 4.1-1 and II, we 
substitute n' ext = —h' for x' 6 <917, (4.22) is reformulated as 

v ex t(x)( 1 - 7u) = + J v ext (x')VG(x,x) ■ fids’ - J G(x, x)X7'v ex t(x) • h'dS' (4.24) 
an an 

This expression is analogous to (4.19), the integrals signs are flipped and y n is substituted by 

7u- 

We now sum (4.19) with (4.24), since the derivative is continuous, (4.4), the second inte- 
grals cancel, reformulating with Vi n t{x) on the LHS we have 


hnt(x) = / rtj nt (x / )V'G(x, x') • h'dS' - 7 r (x)u int (x) 

Jan 

(4.25) 

= [ rtmt(x / )V / G(x, x) • h'dS' + ( 7 u(x) - l)u int (x) 
Jan 

(4.26) 

Finally we can take the last step of our strategy and make 0(x) dependent 
substituting this relation on the potentials split we initially made 

on Uintfax) by 

^( x )mt — u int (x) V vn j (x j 

(4.27) 

= / u int {x')V'G(x , x') • h'dS' + 7 u (x)'Umt(x) 

Jan 

(4.28) 


While initially in order to solve the Poisson problem for cf> with asymptotic boundary condi- 
tions with FEM we needed to bound the surroundings and mesh f7 U I7 C , now with equation 
(4.28) we can solve Poisson’s equation for u meshing only the system 17, from these we have 
the values of Ui n t that we use to compute the boundary values of fa n t ■ We no longer need 
to mesh the surroundings I7 C , because the original Poisson’s problem for (j) has now Dirichlet 
boundary conditions, this greatly reduces the computation cost [30]. 

4.3 Discreet relation between (f) and u 

To obtain the matrix form of (4.28) we can discretize the functions <^>(x) and u(x) using shape 
functions (3.27) and (3.28) from previous Chapter we get 

ne—1 3 ne—1 3 « ne—1 3 

Y Y = ~YY U » / Ar a( x/ )V'G'(x, x') • n’dS' + 7u(x) Y Y 

e=0 a=0 e=0 a=0 e=0 «=0 

(4.29) 

where we dropped the ( int ) notation and assume the mapping is between the internal values of 
4>(x) and u(x) evaluated at x belonging on the surface of the system 17. Choosing a particular 
Xj the LHS becomes fa since only one of the shape functions incident on vertex i is nonzero, 
the others will have a domain that will not contain Xj as nonzero ( cf end of section 3.3). 

Also, on the RHS the terms of double sum that are integrations along a shared tetrahedra 
facets will cancel because h’ will point in opposite directions, so only integrations on the 
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triangles on the surface will be nonzero. Defining an external facet of a tetrahedral element 
close to the surface, d£l c , figure 3.1, as T e ext , (4.29) is reformulated as 

( k = ~ ^2 u a [ ^V'a(x / )V / G(x i ,x / ) • n'dS' + 7u(x i )'u i (4.30) 

{{e,a)\n(e,a)£dfl} Dxt 

This relation can be represented in matrix form as 0 = Gu, where both vectors have as many 
entries as there are vertices within D. A given i will identify a row of G, and all the integrals 
that are coefficients of a given uj will be a part of the entry Gy. If i = j then in addition we 
have 7 u(xj): 

G i,j = - I ^(x')V / G(x i ,x / )-n , d5 / + 7u (x l )0 >j (4.31) 

{(e,a)\n(e,a)=j} xt 

Its our goal now to compute G. 


4.4 Numerical calculation of the integrals and Gy 

Given the elements are tetrahedra, the domain of integration, T^ xt in (4.31) is a collection of 
triangles with different orientations, n! . We can build a single triangle in d = 2 and map it 
to each one of these triangles embedded in d = 3 as in figure 4.3. 

Choosing N 2 points with coordinates (£,7) T within this triangle, each one is associated to 
a given x' G ^% xt - The integrand of (4.31) is then evaluated at these locations and multiplied 
by a weight, the Gaussian weight, and summed, giving us an approximation to the integral. 
Since Gaussian weights are defined for d = 1, further transformations are necessary. The 
figure below shows the mappings we will introduce. After relating a triangle embedded in 
d = 3 with a reference triangle. We will relate the latter with a centered square as in [32], 



Figure 4.3: Sequence of coordinate transformations used to compute the integral Gy. 


4.4.1 Mapping 

Mapping a point from the reference triangle to ( x s , y s , z s ) T G Tl xt can be achieved using shape 
functions, Ai^(x), to represent this location . These functions can be obtained from IV® (x) if 
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n(s, a) = n(e, a) = i and x s G T^ xt 


x s =XqNq(x) + x*Nl(x) + x 2 N$(x.) 

y s =y s 0 NS(x) + yfiVf(x) + y|iV|(x) (4.32) 

z * =ZqNq(x.) + zfiVf(x) + z 2 N 2 (x) 

These three equations are an identity operation, if we give x belonging to an element s as an 
input then the output is still x. Not all three shape functions are independent, noting that 
for every x G fl e we have Nq (x) + iVf(x) + N 2 (x) + /V| (x) = 1 then in particular its also true 
when x G T^ xt . In this case one of the four shape functions will be zero, hence 


iV 0 e (x) + 7Vf(x) + IVf (x) = 1 44 
fl|(x) = l-JV 0 e (x)-IVf(x )44 

^2 ( X ) = 1 — Nq (x) — IV® (x) 


Substitution of (4.33) into (4.32) gives us 


X s =(xq - x 2 )Nq (x) + {x\ - x 2 )iVf(x) + 

V s =(2/o - 2/I)^o ( x ) + (2 A ~ 2/i)^f(x) + 2 /£ 

= (Zq — 4|)/Vq(x) + ( Z 1 - ( x ) + 

Using now the shape functions instead as parameters we define 


Ni 


£ if a = 0 

< y if a = 1 

1 — £ — r) if a = 2 


(4.33) 


(4.34) 


(4.35) 


The reference triangle in figure 4.3 is bounded by y < 1 — £ as shown by noting that 0 < 
Nq 1 2 < 1 implies 


0 < 1 - JVq(x) - JVf(x) < 1 44 

/Vf (x) < 1 - Nq{x) 44 

17 < 1 — £ (4.36) 


This triangle is mapped by the system (4.34), in matrix form we have 


( x s\ 


V s 

V) 

\ 


Xn 


O/Q ^2 ^ 


2/2 2/6 - 2/2 2/l 

y S 
'0 


_ -y 5 



(4.37) 


While x s belongs to a triangular element embedded in d = 3 with vertices Xq, xf and x|, 
(^,r/) T belong to a triangular element defined by y < 1 — £ with 0 < £ < 1 as obtained in 
(4.36). To each point in the last one is mapped into one of the first. When £ = 1 and rj = 0 
then this point inside the reference triangle is associated with Xq, that is the vertex where 
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Nq = 1. In analogous way £ = 0 and i] = 1 will be associated to and when both are zero 
to x|, where iV| = 1 — £ — r]. 


If we know how to map points we can map variations. When £ varies while r] is constant 
results in (dx 3 )^. If instead we vary t] with £ fixed we have (dx s )^. Both are given by: 



( *0 - X 2 \ 


( x\ - x s 2 \ 

{dx s ) v = 

2/o -2/1 

l ^ - 4 J 

d£ and (dx s )^ = 

2/i — 2/1 
\ Z 1 — Z 2 ) 


(4.38) 


With this two vectors we can define a parallelogram. With it, areas are then mapped as 


dS = /dx 5 )^ x (dx s )^| 

= ((Vo4 ~ 4vi) 2 + ( x o4 - 4 x i) 2 + ( x oVi ~ yfci) 2 ) 1/2 d£dri 
= Jd^dr/ (4.39) 


Mapping from a reference triangle to a reference square is represented in figure 4.3. As 
described in [32] first transform the triangle into a square using 


£ = u 

T] = (1 — U)V 


(4.40) 


where 0 < u, v < 1. Variations are given by 


d£\ = / 1 O' 
drj I l — v 1 — u 


du 
dv , 


(4.41) 


and in analogous way areas are mapped by 


dS = (1 — u)dudv 


(4.42) 


Finally we can relate this square to a centered and rescaled one using the transformation 

1 + s 

(4.43) 


u = 


V = 


2 

1 + 1 


with — 1 < s, t < 1, final transformation in figure 4.3. Variations are then as 

/ du\ /l/2 0 \ ( ds\ 

1 dv J lo 1/2 y l dt ) 


and areas are 


dS = -dsdt 
4 


(4.44) 


(4.45) 
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With these transformations, a generic integral along a reference triangle such as 


I = 



W 


f{£,‘n)di i 


(4.46) 


can be reformulated into 


/ = j J f(u, (1 — u)v)(l — u)dudv 

= J 1 J 1 f( 0-m + t) ) l^i dsdt 


8 


(4.47) 

(4.48) 


A Gaussian quadrature as introduced in [12] can now be used, the approximation is done by 
evaluating f(s,t ) at each one of the TV 2 points within the square, multiplied by the Jacobian 
and Gaussian weight WiWj. 


n— 1 7i—l 


i=0 j = 0 



(1 - Sj)(l ~ti)\ 1 - Si 

4 J ~ W ^ 


(4.49) 


Introducing the transformation (4.37) and (4.39) into the integral in (4.31) as well as replacing 
TV® by the definition (4.35) we arrive at 


i-C 


KK T^) 


o Jo 


47t|T s ^ — Xj| 3 


J s d£dr] 


(4.50) 


where the minus sign comes from the gradient of the Green function. This integral is of the 
form (4.46) and can now be approximated analogously substituting into (4.49) 




-Kh(T s d) 


47 t |( T a '£ - x ,;)| 3 


(4.51) 


with 


£ = 
T] = 


1 + S 
2 

(1 - a)(l + t) 
4 


(4.52) 


plugged into 

The next section will show the algorithm to compute these integrals and the G t j entries 
as in (4.31). 


4.5 Algorithm for numerical integration and computation of 
G 

The function boundary_matrix_G in Appendix H will build matrix entries G i j as described 
by (4.31) and the corresponding indices i and j. A given i is chosen and will identify a row of 
the matrix. We then choose a triangular element from the surface s and by running through 
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each a. 6 {0, 1,2} an integral is approximated with (4.49) using N 2 Gaussian points. These 
calculations for a given i are stored as well as the corresponding indices i and j = n(s, a) in 
three separate lists. With these three, the summations in (4.31) are made by csr_matrix for 
the same G t J entry. 

The function is composed by four steps that we will describe below: 

1. Setup Gaussian locations and weights and the Jacobians 

Initially the locations and weights for the reference triangle are built with the function 
gauss_triangle, where Ngauss is given as an input and specifies the number of points along 
one dimension, the function will build Ngauss**2 points locations (£, 77 ) within the reference 
triangle listed in xieta and the corresponding Gaussian weights. 

Then the Jacobians for the transformations given in (4.39) are computed from the locations 
of vertices on the boundary picked from ptot by the indices in alphas. Finally the indice 
lists Ig, Jg for the G entries start empty. 

epsilon=10** (-15) # Tolerance to detect zero integrals . 

xi_eta,weights=gauss_triangle(Ngauss)#I?e/ triangle Gauss Iocs and weights. 
Jacobians=fJ (alphas , ptot) #Jacobians from ref triangles . 

Ig, Jg,G=[] , [] , [] #Where future matrix entries and 

#indices will be kept. 


2. Pick a vertex from the boundary and an element s and test whether the 
integral 4.50 is zero 

A vertex from the surface with index i and location x_i is chosen. A file is initialized on the 
first iteration k=0, or the current indices and matrix entries are saved to prevent these lists 
from getting too large as the for cycles advance. Since in (4.31), (x'(£, 7 /) — x,) ■ n' gives the 
same value for every x.' 6 T^ xt we choose in particular (xq — x* )-h' giving us isnormal, it will 
be used in computing the integrand f later but we can use its value now to decide whether the 
integral will be zero or not. Given that the system we will use is a box, many integrals will be 
zero. If isnormal is zero then further calculations are skipped by the continue preventing 
unnecessary computation. 

for k,i in enumerate (Igs) : #Pick a vertex from the surface with index i. 

" [Code Block] : Initialize a file or save" 

x_i = ptot[:,i] #Location of vertex i. 

for s in range (alphas . shape [1] ) : #Pick a triangle element. 
ig,jg,kg = alphas[:,s] #The 3 vertices locations for s. 
xOs = ptot [ : , ig] 
xls = ptot [: , jg] 
x2s = ptot [: ,kg] 

n = normals[:,s] #The normal to element s. 

d=x0s-x_i #The distance from x_i to element s. 

isnormal=dot(n,d) #If isnormal=0 , skip further calculations . 
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if abs (isnormal) < epsilon: 
continue 


3. Compute the integrals 


If it is not zero then the transformation matrix in (4.37) is computed from the element s 
vertices xOs, xls and x2s. The integrals for each a are computed by summing the integrand 
f for xi and eta along the reference triangle which were computed in step 1. The global 
indices i and j are then saved after being fixed by the block structure of the Jacobian, 
i*5,j*5 gives the 5x5 blocks locations in the Jacobian 5 nq x 5 nq, further addition of +3 
along the columns and +4 along the rows to give the final G entry location. 


#Build the coordinate transformation matrix T 
#from ref triangle to triangle element s. 


xO,yO,zO = xOs 
xl,yl,zl = xls 
x2,y2,z2 = x2s 
T = array ( [[x2,x0-x2,xl-x2] , 

[y2 , y0-y2 , y l-y2] , 

[z2,z0-z2,zl-z2] ] ) 

Js = Jacobians [s] #Pick the corresponding Jacobian. 

for a in range(3): #For a given Ns_alpha compute the integral, I. 

1=0 

for w,(xi,eta) in zip (weights, xi_eta) : 
x = dot (T, array ( [l,xi, eta] )) 
x = x-x_i 

x3 = 4*pi*norm(x) **3 

f = Ns_a(a,xi,eta)*isnormal*Js/x3 

I+=w*f 

G append(-I) #(-l) from derivative of Green fc. 

j = alphas [a] [s] #Append Gif matrix indices fixed 

Ig append(i*5+3) #by the Jacobian block structure. 

Jg append( j *5+4) 


4. Compute the diagonal terms 7 u 

Finally 7u terms are computed for all boundary vertices i. The Kronecker Delta term in 
(4.31) implies that we add these term only to diagonal 5*5 blocks which are then fixed by 
+3 and +4 as in step 3. 

for i in Igs: 

x_i = ptot [ : , i] 

G . append (f gamma (x_i ,pbox) ) 

Ig append(i*5+3) 

Jg append(i*5+4) 

" [Block of code] : Save diagonal entries" 
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The final sum of all contributions for each G % j is made by csr_matrix which will sum all 
those entries in G list that have the same coordinates in IgG and JgG. 


4.6 Restatement of Poisson’s problem 

The original statements that determine (p are now replaced by 


1. Using FEM, find u such that 


V 2 u = V • m if x G 0 


< u = 0 


dujnt 
< dn 


= m • h 


if x € P c 
if x e <9P 


(4.53) 


2. Map u int into (p int using 


(pint — f 1 


(4.54) 


3. Using FEM, find (p such that 


For the <f> residual we have 


\7 2 (p = V • m if x G fi 

<Kx) = (pint if X e <9P 


= / V • ™K(x)d 3 x + J V<j> • S7N*(x)d 3 x 


G2 e 




(4.55) 


(4.56) 


where we do not have the boundary integral since we imposed Dirishlet Boundary conditions 
on <9P [21, p. 92, 108]. Instead of (3.36) we only have now the first two integrals 


s 


e 

(f),a 


src.e 
= s <fi,a 


+ S 


stf,e 

<f>,a 


(4.57) 


The first is given by (3.38) and the second is (3.41) with the IMR substituted. This expression 
is only valid for vertices inside P, that is n(e, a) £ P \ dP. 

The two Jacobians and are now only given by the first terms of (3.57) and (3.56), 
since the second ones are related to the boundary integral. 


dtfy 

n e OS ^a 

dm* 


_ 1 K « 

2 K °‘,P 

P e 

-~S VN ‘ 


(4.58) 

(4.59) 


From the mapping (pint = Guj n ^ for the boundary vertices we can build the condition 0 = 
Guj n t — (pi n t that also has to be verified just like the LLG and the two Poisson’s equations. 
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To evaluate how the boundary values of Mint and 4> tn t satisfy this condition we define a new 
residual just for the boundary vertices 


s 


a 

<t> 


G Hint 4*int 


(4.60) 


each component of these vector is given by 


with IMR: s^ t 


£ 

jGdfl 


Uj Uq.j 
2 


<!>i + 4>0,i 

2 


for i G dQ 


(4.61) 


where we introduced IMR. This is the residual at boundary vertices instead of (4.57), hence 
instead of having the entry (4.58) and (4.59) we have the Jacobians 


ds l,i 

du k 

ds ji 

d<t>k 


G* 


i.k 


2 

1 

~2 


(4.62) 

(4.63) 


which are easily obtained since (4.60) is linear and where i,k € d£l. As in the previous 
Chapter we define the weak form for the residual of u similar to (3.32) 


s* = f V-m N*{x)d 3 x+ [ Vrt • ViV®(x)d 3 a; — f N^Vu-hdS (4.64) 

4n e Jn e Jdti e 

which is split in three parts, the source term is the same as (3.40), the stiffness term is also 
the same as in (3.43) and the surface integral the same as in (3.55) which was obtained 
from the discontinuity of the derivatives in in (3.49). The difference now is that the second 
integral in (3.48) does not exist as we are on the internal boundary and in the first integral 
is substituted the Neumann boundary condition in (4.53) therefore yielding (3.49) as before. 
As a consequence the Jacobians for s® a are given by the two terms in (3.56) and only the 
first term in (3.57). 

Now we have all the residuals necessary to build the system of equations for NR, that is 
the task for next Chapter. 
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Chapter 5 


Residual Vector and Jacobian 
Matrix Assembly 


In Chapter 3 we derived expressions for the linear approximations of the residual of LLG 
and Poisson’s equation r® and s^ a . In Chapter 4 we modified Poisson’s problem to use only 
information from the system by introducing a new potential and a new residual s® a . Each 
one of these three residuals can now be specified for each vertex i yielding linear equations. 

In this Chapter we will show how to organize this set of equations and develop a Python 
code to build it. Of all contributions, we will focus on the two Poisson’s Equations with BEM 
coupling and the exchange term as they will encompass all the main algorithm strategies, the 
other terms are given in Appendix. 

We then introduce our implementation of Newton-Raphson method and briefly refer to 
the algorithms of GMRES with the ILU and Algebraic Multigrid preconditioners to solve the 
linear system of equations. 


5.1 A larger Jacobian, J 


The configuration for every vertex in the system 0 is specified by the magnetization m.; 
and two potentials <f>i and u t which can be assembled into a vector with five entries x = 
(m x , m y , rn z , 0, u)f . To measure how close a local configuration composed by the vertex i and 
its connected neighbours is to the exact solutions of LLG equation and Poisson’s equations 
for both potentials, we define at i a five entry residual vector yy = (r x , r y , r z , s 0 , s u )J. 

If we have a particular configuration for i and his j neighbours, if we change the config- 
uration for each j, how does the residual for i change? Consider for example the first order 
approximation of the discretized LLG residual. When only the magnetization nx,- varies from 
a particular configuration, m j tP , the variation of the residual at i is obtained from the 3x3 
Jacobian matrices of first derivatives and the following linear approximation 



/ \ 


j dr x 

^x,p 1 


1 dm x 

r y,p 

+ 

dry 

dm x 

\ r z,p) 


l dr. 

i 

\ dm x 


dr x 

din v 

dry 

dirty 

dr , 

dirty 





( A m x 
A m y 
\A m z 


+•••+ 


i 


/ dr x 
dm x 
dry 
dm x 
l dr z 
\ dm x 


dr x 

dm v 

dry 

drily 

dr z 

dirty 



dntz I 

dr, I 
dm, / j j 



+ . • • 
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These two matrices shown represent the relation between the residual at vertex i and a 
variation in magnetization at i itself and at a neighbour j, their expressions were obtained 
in previous chapters for each term of LLG residual when two vertices are connected. When 
c t>j varies as well an additional vector term is added to (5.1) given by (f>i- If 

there is no edge in the mesh connecting (i,j), this matrix and vector are made out of zeros. 

An analogous expression for the linear dependence of or s Ut j with magnetization and 
potential at j is given by 


S(j), i 


+ * * ‘ + 


V dm x,j 


d m y,j 


dm z } j 


( A m x 
A m y 
\A m z 


, ds *j A A 
+ w A cpj + 


dcf> 


j 


(5.2) 


except when i G dQ where s ^ j has the form (4.61). 

We can organize the dependence of the four residuals with the variation of a configuration 
of a given j neighbour as 


^ r x \ 

r v 


r y,p 


( 

dr 

dm 

dr 

d<)> 

0 

0 

r z 

= 

^ z iP 

+ 




0 

S(f) 


S(j),p 


dm 

ds ^ 
d<t> 

0 

\s u ) 

i 

J 

i 


ds u 

dm 

0 

ds u 

d<j> 


^Am^ 

A m y 
A m z 
A <j> 

\ Au h 


(5.3) 


This 5x5 matrix defined as J y is associated with a generic edge (i,j) of the mesh, it has 
several blocks. The is a 3 x 3 matrix as in (5.1) that is composed by the time derivative 
term A® a as in (3.64), the triple term contribution given by D® a as (3.76) and the exchange 
contribution E® ^ given in (3.94). The 3x1 column vector is V e a j3 in (3.77). The two 
1x3 row vectors and both defined as Q® ^ in (3.56). And the the two scalars ^ 
and ^r- both derived in (3.57). The variation vector in (5.3) is defined as Axj and multiplies 
the matrix giving how much the residual changes when the particular configuration for the j 
vertex changes as well. 

Since i has several neighbours we have to extend (5.3) to include their Jacobians. Therefore 
in an analogous way as we summed matrix vector multiplications in (5.1) we have 


T‘i — T Ji,oA)Co + • • • + JijAXj + • • • T Ji,n 9 — 1 Ax Uq —\ (5'4;) 

By joining all 5 x 5 blocks in (5.4) into a 5 x 5 n q matrix and all n q variations into a single 
5 n q x 1 column vector this sum can be represented as a dot product of both, which represents 
all possible relations of i with all vertices in the mesh including vertices that are not connected. 
This situation is analogous to what we have did to obtain (3.45) but with a much larger matrix 
and vector. 

We can go further now and establish the same relation (5.4) for all i vertices generating 
n q matrices with shape 5 x 5 n q . Observe now that each one of these matrices is in fact a 
5 x 5 n q row block of a much larger Jacobian matrix J which results from assembling all these 
into a final 5 n q x 5 n q Jacobian. By setting the LHS of equation (5.4) equal to zero for all i 
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and rearranging we arrive at 



This is the linear system of equations that we have to solve iteratively as we described in 
Section 3.2. Remember that J on the LHS and the RHS residual, r, both depend on mo, (f ) o 
and uq and then are evaluated at a particular guessed configuration, that we have decided to 
start at mo, <f > o and uq. After assembling the J and r we solve for the variations, Ax. We add 
this variation to our guess giving us the next particular configurations we use to recompute J 
and r while maintaining the same mo, (f > o and uq. As we proceed, the residual will on the RHS 
get smaller progressively until a prescribed tolerance. The current particular configuration 
thereby obtained is close to the exact solutions, that is we have an approximation for the 
update of the initial magnetization and potentials that satisfies the LLG equation and both 
Poisson’s equations. Repeating the whole process we construct the history of the system from 
the initial configuration. 

Notice we have to assemble and recompute some of the blocks of this huge sparse matrix, 
J, and vector, r, at every iteration. This has to be done quickly. In the next section we review 
the techniques proposed by [24] for the assembly of Poisson’s stiffness matrix and we propose 
the extension of such techniques for all terms in the LLG equation. 


5.2 Poisson’s Contribution 


To build the Poisson’s stiffness matrix we have to return to the gradients of the basis functions 
and introduce the technique to compute them. Remember from (3.23) that for an element e 
the basis has the form 


^a( x ) — X7vr( c a,0 + c a,l x + c a,2V + c a,3 z ) 


where c e a i is the cofactor of entry (a, i) of the matrix X composed by the spatial locations for 
each vertex in a given element. Its gradient is given by 


C a ,i\ 

= C “- 2 
\Ca, 3 / 


(5.7) 


and is independent of x. 
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To build the stiffness matrix we need to compute these gradients for all elements and that 
in turn requires computing the three cofactors for each one. From (3.20) take for instance 

yi zi 

2/3 z$ 

the minus sign outside the parenthesis is the signal of this cofactor. A simple algorithm to 
compute all gradients would be to pick an element then its four vertices’ global indices on 
column e of the connectivity matrix, whose spatial locations are used as inputs for each one 
of the three cofactor expressions such as (5.8), giving us one gradient for a particular (e,a). 
Repeating the same process for all remaining a we get the four gradients of the four basis 
function of an element. These calculations are then carried on all other elements and saved. 

A severe drawback for an algorithm like this in Python is the use of many for cycles, 
which are very slow [24], In particular the issue is not mainly related with the Poisson’s 
stiffness term for (j) and u since it is built out of constant gradient terms and hence only needs 
to be computed once. The issue is when we have to go through the process of computing and 
assembling Jacobian blocks repeatedly, as occurs for all other terms of LLG which changes as 
magnetization and potentials evolve. 

As [24] points out, computing the gradients and then the K entries in (3.44) can be 
reformulated with vectorial operations such as element-wise multiplications and sums, which 
by themselves are faster and reduce the number of required for cycles. The central idea is to 
vectorize all possible operations. 

In the next section we will introduce these concepts for Poisson’s stiffness matrix from [24] 
and then go further and extend them to develop the code for the J part associated with the 
Poisson’s stiffness matrix for both potentials taking into account the BEM mapping between 
them, we build the three vector residuals associated with (3.36) for both potentials as well, 
that, once summed, constitute one of the parts of the final residual vector r in (5.5). 

5.2.1 From Gradients to Poisson’s Stiffness Matrix 

Our goal is to build the n q x n q matrix entries in (3.44) and then relate them to the J matrix. 
Start noticing that co,i above can be reformulated. Temporarily complicating the expression 
we have 

2/1 z\ ]j2 z 2 , 2/1 zi i/I zi 

On = - + - 

2/1 zi 1/3 ^3 2/3 z$ 1/2 Z 2 

Since determinants are linear row-wise we join the first with the last determinant and the 
middle ones 

2/1 Z! 

2/1 - 2/2 zi - 2:2 

Multiplying the second by (— 1)( — 1) yields 

Vl 2/3 ~ 3 = (2/1 - 2/3) (zi - Z 2 ) - (2/1 - 2/2) (zi - z 3 ) ( 5 . 11 ) 

2/1 - 2/2 zi- z 2 
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Co,l = 



co,i = 


(5.9) 





Unlike the original expression (5.8) which had mixed components, the latter is a product of 
differences of the same component. This is important because it unlocks the possibility to 
use matrix element-wise operations provided we have the right data structures. Hence, we 
maintain the connectivity matrix as an 4 x n e array and we restructure the array of vertices 
locations from the shape n q x 3 to 3 x n q . 

Now, instead of computing all cofactors for a given element and proceeding with the 
next one as suggested in the simple algorithm, what we want is to compute the gradients 
for all vertices associated with a given a as a single vectorial operation. This is achieved as 
follows, we start by defining alphae as a Python array with the connectivity matrix entries, 
noticing that alphae [a] picks the a row, with n e entries, which contains all global indices of 
all a vertices. From the array of locations ptot, all locations of all these vertices are sliced 
as ptot[: , alphae [a]], which is a 3 x n e array where each column gives us the position 
in space of all vertices i = n(e, a = a) for all e. These slices allow us to compute the 
required differences in (5.11), for example the first one, y\ — y 3 , can be computed for all e as 
ptot [1 , alphae [1] ] -ptot [1 , alphae [3] ] where the first index 1 indicates we are choosing 
row one of ptot, the row with all y components. Since we will ultimately require differences 
between all components we do better computing differences between vectors and later pick 
the right components. The following difference 


( x£ c\ 

Va 

W 


4\ 

y e A 

W 


(5.12) 


between the position of a vertices and the position of /3 vertices can be computed across all 
elements using 

D_ab = ptot [:, alphae [a] ] -ptot [:, alphae [b] ] 


Since (5.11) is composed of differences between components it can be computed for all 
elements using the arrays D_12 and D_13. We pick the second row of D_13, the one with the 
differences yf — y\ for all elements, and the third row of D_12, which has the z component 
differences z\ — zf- Both are then multiplied element-wise as D_12[l, :]*D_13[2, :] yielding 
(yf — — zf) for all e, which is then summed to the second term in (5.11) as 

(D_13 [1 , : ] *D_12 [2 , : ] -D_12 [1 , :]*D_13[2, :])*C 


where the constant C=l/ (6*signedvol) has all determinants of X for all elements element- 
wise multiplies every entry of the array in parenthesis. 

The other two cofactor components for V/Vq are the second and third entries of 


VNo e (x) 


1 

det(X) 


( (yi ~ 2/3X21 ~ 22) - (2/1 - 2/2X21 - 23) 
(Xl - X 2 ){z\ - Z 3 ) - (Xl - X 3 )(z\ - Z2) 

\(aq - x 3 )(y 1 - 2 / 2 ) - (xi - x 2 )(yi - 2 / 3 ) 


(5.13) 


We observe that to assemble all gradients for a = 0 can be done extending the element-wise 
operations we used for the first entry to the other two. Defining now a 3 x 4n e array of zeros 
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Grads, we can fill the first 3 x n e block with the gradients ViVg(x) as follows 
Grads=zeros( [3,4*ne] ) #Initial G. 

Grads [: , :ne] = array( [ (D_13 [1 , :]*D_12[2, :]-D_12[l, :]*D_13[2, :])*C, 

(D_12 [0 , : ] *D_13 [2 , : ] ~D_13 [0 , :]*D_12[2, :])*C, 

(D_13 [0 , :]*D_12[1, :]-D_12[0, :]*D_13[1, :])*C]) 

The remaining three block with shape 3 x n e will receive the ViVf(x), ViV|(x) and V7V|(x) 
calculations 


Grads = 


a = 0 

a = 1 

a = 2 

a = 3 

[ :ne] 

[ne : 2*ne] 

[2*ne : 3*ne] 

[3*ne : ] 


(5.14) 


has the Algorithm A. 4 shows [24]. 

Once we have all entries for the matrix Grads, its structure allows fast computing of the 
stiffness matrix entries. Remember that each (i, j) entry of K matrix is given by 


Kij = ^ VA 1 ‘ VIV^ 6 (5.15) 

{(e,a,/3)\i=n(e,a)/\j=n(e,/3)} 

We observe is a summation of dot products, p, mapped by the connectivity matrix into 
the same All entries of K then will require all a calculations, we start by organizing 

them as consecutive blocks of size n e in the following Kg array 


' (0,0) 

(0,1) 

(0,2) 


{a, f3) 


(3,3) 

[:ne] 

[ne : 2ne] 

[2*ne : 3*ne] 


[c*ne : (c+1) *ne] 


[15*ne : ] 


— — X lU/fcg 

(5.16) 

where c stands for the index in {0, 1, ..., 15} associated to a combination (a,/3). Before han- 
dling the details of how to compute each entry in Kg we will focus now on its block structure 
to build the indices (i,j) for each entry. Defining Ig and Jg as two arrays of zeros with the 
same shape as Kg and the same 16 row block structure with length n e 


Ig 


0 0 0 


a 


3 


lx 16n e 


(5.17) 


Jg 


0 12 



3 


lx 16n e 


(5.18) 


for each Kg entry we assign for both the respective global index i = n(e,a ) and j = n(e,j3). 
This is done with two simple arrays ii and jj [24], associated to a and (3 sequences in the 
n e blocks of Kg 


ii 


0,0, 0,0, 1,1, 1,1, 2, 2, 2, 2, 3, 3, 3, 3 


J 1x16 


(5.19) 


JJ = 


0 , 1 , 2 , 3 , 0 , 1 , 2 , 3 , 0 , 1 , 2 , 3 , 0 , 1 , 2, 3 


J 1x16 


(5.20) 


We cycle through all (a,/3) combinations by ranging c from 0 to 15 and picking a = ii[c] 
and j3 = j j [c] , then we fill each n e block initially with zeros in Ig and Jg respectively with 
the n e global indices associated to a and (3, that is with alphae [ii [c] ] and alphae [j j [c] ] . 
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Once we have the three arrays, Kg, Ig and Jg, we give them as inputs to the csr_matrix 
Python function that builds the matrix and entries and if two or more entries have the same 
combination of (i,j), they are summed and assigned into that entry in the matrix. This 
function yielding an array in csr format and shape n q x n q . 

This array could be used in solving the Poisson’s equation alone. However our goal is 
to make it a part of the larger Jacobian matrix. Observe that each J lj block in J has the 
structure given in (5.3). The entry K t j is equal to the two Poisson entries in the J t .j block 
in J. 

Where for the moment we ignore the fact that the boundary vertices have the residual 
given by (4.60). Observe that in J the global indices 5 i and 5 j locate the upperleft entry 
of Ji.j 5x5 block. To locate the Poisson’s < ft entry in J we need to add a local increment, 
(5* +3, 5 j +3), for potential u we have (5i+4, 5 j +4). This observation amounts to multiplying 
Ig and Jg by 5 and then adding a local fixing +3. Therefore, to build the indices arrays that 
map the Kg 1 x 16n e entries into the Poisson entries in J we proceed as 

Ig = zeros( [16*ne,] ,dtype=' int 1 )#Initial Ig and Jg 

Jg = zeros ( [16*ne ,] ,dtype= 1 int 1 ) 

for c in range (4*4) : #Choose a combination (alpha, beta) 

Ig [c*ne : (c+1) *ne] =alphae [ii [c] , : ] *5+3 #(+ 4 ) if u 
Jg [c*ne : (c+1) *ne] =alphae [j j [c] , : ] *5+3 #(+ 4 ) if u 


The details to finally handle each entry in each n e block of Kg we will use a particular 
example. Choose the first two columns of Grads, each one is associated with element e = 0 
for vertices a = 0 and /3 = 1, the dot product is 


VJV£ =0 • ViV° =1 = 


1 


(6fi e 


,0 


„o 


e\2 V 0,a: ^0,y 



— c 0,x c l,x + c 0,x c l,y + c 0,z c l,z (5-21) 


A careful observation shows this operation is the same as multiplying both column vectors 
element-wise and then sum the vector column-wise. The structure of Grads is the one that 
allows a generalization, we element-wise multiply two 3 x n e blocks associated with a and (3 
resulting still in a 3 x n e array, that once summed column-wise yields a row vector 1 x n e 
where each one of the entries is ViV® • ViV|fl e . 

Hence the first five (a, {3) combinations that will be used in computing the stiffness matrix 
as in [24] are given by 


Kg = zeros( [ne*16,] ) 

#(alpha, beta) 

Kg [0 :ne] 

= sum(G_0*G_0) *vol 

#(0,0) 

Kg [ne : 2*ne] 

= sum(G_0*G_l) *vol 

#(0,1) 

Kg [2*ne : 3*ne] 

= sum(G_0*G_2) *vol 

#(0,2) 

Kg [3*ne : 4*ne] 

= sum(G_0*G_3) *vol 

#(0,3) 

Kg [4*ne : 5*ne] 

= Kg[ne:2*ne] 

#(1,0) 
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where we observe that since K is symmetric the (0, 1) block is the same as the (1,0). All 
other ct, j3 combinations are in Algorithm B.4. 

Now that we have the necessary information to build both and for the (i, j) 
blocks of J, using the Kg array together with indices Ig and Jg we have 

Kgphi = csr_matrix( (Kg, (Ig, Jg) , shape=(5*nq, 5*nq) ) 


as a 5 n q x 5 n q sparse matrix where only the entries within the 5x5 blocks associated with 
connected vertices are nonzero. In a similar manner, the process of assembling the stiffness 
matrix for the potential u is given by 

Kgu = csr_matrix( (Kg, (Ig+1, Jg+1) ,shape=(5*nq,5*nq)) 


since the entries are the same but one entry further down the diagonal of each 5x5 block. 

However, remember that for the vertices that belong to we have the potential cj)i 
determined by all boundary values of the potential Uj and not an equation that related it 
with the neighbours represented by the coefficients at each row for all i within H. We have 
then to replace all entries in Kgphi at every row associated with i E <912 by (4.63) and (4.62). 
In practice we define first the Igs as an array of global indices of vertices at the boundary, 
second we assemble the array G introduced in section 4.5 that maps u into </>, then we introduce 
the Jacobians of (4.62) as 

Kgphi [Igs*5+3 , Igs*5+3] =-l 
Kgphi+=G 


where we postpone the introduction of the \ factor from IMR. By summing Kgphi and Kgu 
we will have a 5 n q x 5 n q csr matrix that can later be added to other assembled Jacobian 
contributions. 

An important observation from the two previous contributions is that they do not depend 
on the current configuration x, meaning it has only to be done once. In fact there are other 
contributions that are also constants as well. The time derivative term A in Algorithm F.l and 
both Q for u and </> are also constant and described in Algorithm E and can be added to the 
previous two Poisson’s terms. We will call it Jfix and each one of these terms is computed 
and assemble by the function J_f ix_assembly in Algorithm G.l, note that, apart from the 
time derivative, the \ factor from IMR is present in every Jacobian term, hence only after we 
sum all these terms 

Jfix = (Kg_phi + Kg_u + Q_phi + Q_u + Q_u_BC)*0.5 + A 


we multiply by this factor and only then we add A. Where we also define 
KgphiuG = Kg_phi + Kg_u 

with both -I and G included, it will be useful to build the residuals next. 
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5.2.2 Assembling the two Poisson’s Residuals 


Each one of the two Poisson’s residuals has three integrals contributions in (3.36). In Chapter 
3 we derived their discretizations, giving us three 5 n q x 1 vectors for which we now develop 
the code for a single residual vector for both potentials. 

The first of the three terms is the source term, it requires first to define a new structure 
for the magnetization. The column vector x has each magnetization is at every x [5*i : 5*i+3] 
slice, from these we reorganize them into a 3 x 4 n e array malpha, where each 3 x n e block is 
associated with an a and each column is . 


malpha 


a = 0 a = 1 


a = 2 


a = 3 


3x4n e 


(5.22) 


obtained in A.l. From (3.40) we will require the divergence of the magnetization (3.37) for 
each element, it is given by the sum of the dot product between the magnetization, m® , with 
the gradient of the basis function, ViV®, for the four a, each of these can then be repeated for 
all e. We can accelerate this procedure by changing what we assume as being fixed. Instead of 
fixing e and cycling through a we are better fixing a and vectorize the calculations for all e. 
Therefore we multiply element-wise malphaf : , a*ne : (a+1) *ne] *Grads [ : , a*ne : (a+1) *ne] 
for a given a = a giving us a 3 x n e array, if now we summed this array column-wise it would 
yield V • (m e a Nff) for every e, but if we first sum the former for all other a, then we will only 
need to make a single sum 

divm = sum(malphaO*GradsO+malphal*Gradsl+ 
malpha2*Grads2+malpha3*Grads3) 


giving us a 1 x n e array with the divergences for every element, (3.37). 

To build the first contribution for the residual we will use this divm array. Since the 
divergence only depends on the element e, the calculation in (3.40) is obtained by summing 
all divergences for all elements that have i as one of its vertices, this is done by fs_src_phi 
function of Algorithm B.2 together with the indices in B.3, giving us an array s_src_phi_u 
with shape 5 n q X 1. 

The corresponding residual source term for the u potential is exactly the same as the one just 
computed. But one entry below in the 5 n q x 1 vector. Therefore to build a residual vector 
for both Poisson’s problems, we can simply slice each fourth entry at every 5i x 1 block and 
map them into the u entries as s_src_phi_u[4: : 5 ] =s_src_phi_u [3 : : 5 ] as we do in #1 in 
the following algorithm, B.6 

#l.The source term. 

s_src_phi_u=f s_src_phi (malpha , vol , Grads , nq , barcode) 

s_src_phi_u [4 :: 5] =s_src_phi_u [3 :: 5] #copy src term of phi to u entries. 

s_src_phi_u [Igs*5+3] =0 #zero the phi boundary entries. 

#2. The boundary u term. 

s_u_BC=f s_u_BC (malphas , alphas , normal_vecs , areas , nq) 

#3 .Compute the residual (KgphiuG does not have 0.5, its in avg x) 

sphiu=s_src_phi_u+KgphiuG . dot (x) -s_u_BC 
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However the array s_src_phi_u still needs further modification, since for every i vertex at 
the boundary we do not have a Poisson’s residual but the residual in (4.60). Hence, for these 
5i x 1 slices of s_src_phi_u we do not have the entries just computed, instead of changing 
fs_src_phi and its indices we simply set them equal to zero as s_src_phi_u [Igs*5+3] =0. 

The residual at these vertices can be joined with the stiffness residual vector for both po- 
tentials and the BEM coupling. This vector term is essentially the KgphiuG matrix computed 
in #2 times the average of configuration x=(x+x0)*0.5, note the 1/2 IMR factor is included 
on the average configuration and not in the Jacobian itself. 

The final integral for the Poisson’s residual is only present for the u potential, it is com- 
puted by the f_s_u_BC function in Algorithm B.5, together with the stiffness terms and 
residual for the boundary vertices, the residual for the two Poisson’s problems the 5 n q x 1 
column csr array sphiu that corresponds to s from Section 4.6. 


5.3 Exchange Matrix Block 


We will focus on the calculation of (3.93) and only then substitute IMR. It can be split into 
two parts. We will develop first 
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E 


m /3 



(5.23) 


This is done by the function Jex in C.l. We will combine malpha structure with the Poisson’s 
entries Kg in (5.16). Consider a particular a, then (5.23) is a single sum along the /3 € 
{0, 1, 2, 3}. It can be extended to every e by 


1 

4 


3 


E 


m 

m 


o 


0 

y,P 


m 


o 

z,(3 


m 

m 

m 


n e —l 
x,/3 
n e — 1 

v,P 

n e — 1 

z,P 


zv° 


Ty-n e — \ 


(5.24) 


with both blocks having n e columns. The first is a slice of (5.22) given by 
malpha [: ,b*ne: (b+l)*ne] and the second, the slice of Kg associated with (a,/3). We mul- 
tiply each column of the first by the corresponding entry of the second, giving us a 3 X n e 
array which is summed element-wise with the others for f3 = 1,2,3. The calculation (5.24) 
have now to be extended for all a. Starting with a 3 X 4n e array BO with only zero entries 


BO = zeros ( [3,4*ne] ) 
for c in range(16): 

a=ii[c] #alpha 

b=j j [c] #beta 

BO [ : , a*ne : (a+1) *ne] +=malpha[ : ,b*ne : (b+1) *ne] *Kg [c*ne : (c+1) *ne] 


we compute each summand in (5.24) and add it to the right block in BO which is guaranteed 
by the sequences in ii and jj introduced previously for Poisson blocks. The final result is 
a BO array with four consecutive blocks, each one with shape 3 X n e and representing the 
calculation given in (5.24) for a given a. 
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K e 

The second term is given by mf ot ^ ,7 . For a given e, each combination c in {0, 1, 15} 
determines the „ while mf ot is constant. Defining a mtot as 3 x n e array with columns as 
(3.84) we compute it as Algorithm A. 3. Using just the sequence ii this single multiplication 
is extended to all e as Kg[c*ne: (c+1) *ne] *mtot which is then added to BO sliced at the 
corresponding a 

Bvec = zeros ( [3, 16*ne] ) 
for c in range(16): 

a=ii[c] #alpha. 

Bvec [ : , c*ne : (c+1) *ne] =(B0 [ : , a*ne : (a+1) *ne] - 

Kg [c*ne : (c+1) *ne] *mtot) /4 


the 1/4 factor in (3.93) is only now introduced since the BO term will be used in the exchange 
residual calculation without it. The final array will be structured as 


Bvec = 


(0,0) (0,1) (0,2) ... (a,/3) ... (3,3) 


(5.25) 


Remember that building the Poisson’s matrix required our data aligned in a single row, so 
that two arrays of indices could map all entries into the final Jacobian matrix. Our array 
Bvec has three rows and also notice we now have to map to a 3 x 3 matrix block embedded 
in a 5 x 5 block with indices (i,j) instead of a single Poisson entry. We will divide the 3x3 
block into its lower triangular part and upper triangular, and then build the indices for both, 
but first Bvec has to be transformed into an array with shape 1 x (3 * 16n e ). Notice in (3.75) 
the entries for the lower triangular part have the y component negative, hence 

B=zeros( [16*ne*3*2,] ) 

#Lower triangular. Flatten and app (-1) to y components 
B[:16*ne] = Bvec[0,:] #x 

B [16*ne : 2*16*ne] = -Bvec[l,:] #y 

B [2*16*ne : 3*16*ne] = Bvec[2,:] #z 

#Upper triangular has oposite signal. 

B[3*16*ne:] = -B[:3*16*ne] 


Since its a skew-symmetric matrix block, the upper triangular part is the minus the lower 
one, we then copy all previous entries, B [ : 3*16*ne] and multiply them by -1 and add to 
B[3*16*ne:] . The result is a single row array B divided into two blocks, and each one has 
three sub-blocks of length 16n e , associated with the x, y and z components and in turn for a 
given component we have a (a, j3) structure with 16n e entries, that is one of the rows of Bvec. 
From another point of view, the idea is to flatten Bvec x, y and z rows twice, and ’’glue” both 
arrays, the difference between them are the signals, while on the first half the y components 
have negative sign and x and z have no signs, on the other its the opposite. 

What remains to be done from this structure and the structure of matrix J is to build 
the corresponding indices for B array. We start bottom up and observe that each row of Bvec 
has the same (ct, f3) structure as Kg therefore we build the indices from alphae rows similar 
to what we did for Poisson’s term (p. 59). 
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#l.For x row of Bvec build the sequence of 
# (alpha (i) , gamma (j)) indices 

ii=outer ( [0 , 1,2,3] , [1,1, 1,1]) f lattenO 
j j=outer ( [1,1, 1,1] , [0, 1,2,3]) f lattenO 
#2. For a row of Bvec build the (16) block indices 
Ig = zeros ( [16*ne,] ,dtype=' int 1 ) 

Jg = zeros ( [16*ne ,] ,dtype= ' int 1 ) 
for c in range(16): 

Ig [c*ne : (c+1) *ne] =alphae [ii [c] , : ] 

Jg [c*ne : (c+1) *ne] =alphae [j j [c] , : ] 


Then we copy both Ig and Jg three times, one for each row of Bvec. We multiply it by 5 to 
yield the locations of the upper left entry of the (i, j) block of 5 x 5, and then add +1, or +2 
or +3 local fixing depending on whether we are building the indices for the lower triangular 
or upper triangular. Explicitly we add them as follows 

IgEx=zeros( [3*16*ne*2,] ,dtype=' int ' ) 

JgEx=zeros( [3*16*ne*2,] ,dtype=' int ' ) 

To the 16 block locations in Jacobian, 

#add the local x,y,z increments 
Igxyzf ix=array ([2,2,1]) 

Jgxyzf ix=array ( [1,0,0] ) 
for d in range (3): 

#lower trianglar 

IgEx [d*ne*16 : (d+1) *ne*16] =Ig*5+Igxyzf ix [d] 

JgEx [d*ne*16 : (d+1) *ne*16] =Jg*5+Jgxyzf ix [d] 

#upper triangular 

IgEx [(3+d) *ne*16 : (3+d+l) *ne*16] =Ig*5+ Jgxyzf ix [d] 

JgEx [(3+d) *ne*16 : (3+d+l) *ne* 16] =Jg*5+Igxyzf ix [d] 


Where Igxyzf ix[d] and Igxyzf ix[d] entries are with respect to the lower triangular part, 
for example d=l corresponds to a y entry, therefore from (3.75) we need to add +2 down the 
i direction while j direction does not require fixing, hence +0. For the upper triangular part 
we switch them. The final result is three arrays B, IgEx and JgEx which serve as input to 
csr_matrix. The J matrix contributions from exchange term is then added to the Poisson’s 
matrix, notice both will be 5 n q X 5 n q matrices and both are pieces of the final Jacobian. 

In order to include IMR, this function needs the input variable, malpha, to be the average 
between malphaO associated with all 0 and the current guess m^ p , malpha. 

5.3.1 Exchange Residual 

Observe in (3.85) that the term in parenthesis is the same as the first term in (3.93) which 
was already computed in B0, for a given a and e we compute the cross product of m® oi with 
m e pK^ p. We can extend this operation for all elements, using B0, since it is already 
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organized in four blocks of 3 x n e one for each a. The total magnetization for each element 
constitutes the 3 x n e array mtot. What we do now is to use these data structures and 
introduce a vectorized version of the cross product, this is done by noticing that 


"a x b] i — (ij bk 0 k bj 


(5.26) 


implies we can compute the cross product by summing two element-wise multiplications be- 
tween the vectors a and b with the entries permuted as follows 



a\ 
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(5.27) 


The next lines of code will compute these cross product for a E {0, 1, 2, 3} for all e 

#1. Cross is computed for each vertex alpha 
rvec=zeros ( [3 , 4*ne] ) 
for a in range (4): 

x = mtot [ [1 , 2,0] , : ] *B0 [ [2, 0, 1] , a*ne : (a+1) *ne] #Permute rows in mtot 
y = mtot [ [2, 0, 1] , : ] *B0 [ [1 , 2,0] , a*ne : (a+1) *ne] #and B0 and e-wise mult. 
rvec [ : ,a*ne : (a+1) *ne] = -(x-y)/4 #-(cross product )/4- 


Where for example mtot [ [1 , 2 , 0] , : ] is a 3 X n e array obtained from mtot by flipping the first 
and second rows and then the new second and third. Similarly B0 [ [2 , 0 , 1] , a*ne : (a+1) *ne] 
block has its rows flipped as [ 62 , bo, bl] T in (5.27). The output is an array with shape 3 x 4n e 
where each column is given by (3.86) for a given a. 

In order to map these vectors into the final residual vector we have to flatten these three 
rows in similar manner as we did previously with the exchange Jacobian. 


#2.Flatning 



r_EX=zeros ( [3*4*ne 

,]) 


r_EX[:4*ne] 

rvec 

[0] 

r_EX [4*ne : 8*ne] = 

rvec 

[1] 

r_EX [8*ne : ] 

rvec 

[2] 


The indices that will map each entry of rEX vector into the final 5 n q x 1 vector are built 
by observing the structure of one of the rows of rvec, it has the same a block structure as 
malpha. We start by building an array of indices Ig for this structure as 

#1 . Define indices for the x component . 

Ig=zeros ( [4*ne , ] , dtype= 1 int ' ) 
for a in range (4): 

Ig [a*ne : (a+1) *ne] =alphae [a] 


Noticing that r_Ex has a structure given by three consecutive flattened malpha structures 
giving us a 1 x (3 * 4 * n e ) array. We modify Ig to give us the final indices by multiplying by 
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5 for the locations of the first entry in each 5x1 block of the final residual, onto which we 
then add +0 if its an x entry, +1 if its y and +2 if its a 2 entry depending on which one of the 
the three blocks of indices we are building. 

#2. Build indices for each one of the three rows of 
#rvec in rEx function. 

Igxyzf ix=array( [0, 1 ,2] ) #x, y ,z components 
IgEx=zeros ( [3*4*ne , ] ) 
for d in range (3): 

IgEx [d*4*ne : (d+1) *4*ne] =Ig*5+Igxyzf ix [d] 


5.4 Assembling all terms 


The final algorithm will assemble all matrix J and the residual vector r from an initial con- 
figuration xO and a particular one x. 

The J blocks contributions are already available now by the functions previously described 
and all the other in the Appendix. Hence the remaining task is to simply initiate each of these 
functions in sequence, but first we have to build their inputs from x and xO, we will require 
m. e a — m® 0 , and the averages for m® , 0® and mf ot a , which are computed as follows 


avg=zeros ( [3 ,nq] ) 
dx=x-x0 


#1 . Initiate an array with zeros. 

#2. Compute the difference between configs. 

5] #3. Put each component of mag at 

5] #one row of avg, the natural 

5] #structure for fmalpha. 

dmalpha = fmalpha(avg, alphae) 

x=(x+x0)*0.5 #f. Compute the avg (IMR) config and replace x. 


avg [0] = dx [0 

avg[l] = dx[l 
avg [2] = dx [2 


avg [0] = x [0 

avg [1] = x[l 

avg [2] = x [2 


:5] #5. Build an array where 

:5] #each column is the mag 

:5] #at i. 

avgmalpha = fmalpha (avg, alphae) #6. Mag at each (e, alpha) . 
avgmalphas = fmalphaS (avg, alphas) #7. Mag at each (s, alpha). 


avg [0] = happO [0 : 
avg[l] = happO [1: 
avg [2] = happO [2 : 


3] 

3] 

3] 


#8. Do the same for the cte 
# applied field. . . 


avghappalpha = fmalpha (avg, alphae) 

avg[0] = x[3::5] #9... and the potential phi. 

avgphialpha = fphialpha(avg [0] , alphae) 
avgmtot = fmtot (avgmalpha, ne) 


avg [0] = xO [0 

avg[l] = x0[l 
avg [2] = xO [2 


5] 

5] 

5] 


#10. Malpha structure for mO. 
#No average used. 


malphaO = fmalpha(avg, alphae) 
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where we recycle the array avg through the sequence. 

Finally on top of Jf ix that was already loaded from memory we add all other contributions 
as in Algorithm G.2 yielding the J 5 n q x 5 n q matrix and column array r that are an essential 
part for the Newton-Raphson method we revisit on the next and final chapter. 


5.5 Newton-Raphson algorithm 

From sections 4.6 and 3.2 we will build our code for the algorithm. We will define two simple 
functions. The NewtonRaphson function in Algorithm 1.2 starts with the configuration xO 
and the particular configuration x=xO and solves a sequence of linear systems to yield the 
configuration x at t + At. And the Update function in Algorithm 1.1 that simply saves this 
x configuration and restarts the NewtonRaphson function for the new x. Additionally if the 
applied magnetic field varies with time that process would also be here introduced, but we 
will choose to keep it to be constant. 

The first iteration of the while loop below defines the Jacobian and residual of the linear 
system of equation (5.5) evaluated at a particular configuration x equal to xO computed by 

Jandresidual. 

while 1: 

J,r = Jandresidual (x,xO, hap, hapO,deltat,normal_vecs, 

alphas , alphae , Igs , alphac , Grads , areas , 
vol , nq , ne , Jf ix , KgphiuG , barcode) 
if it>it_max or norm(r) <restol : 
return x 

deltax=rILU_AMG_GMRES (J, -array (r) reshape (5*nq, ) , 

array (xO) , args_phi , args_u , args_gmres) 

x+=deltax 


The residual norm is then computed and compared against a given tolerance, restol. Initially 
this step is skipped as the residual is higher than the tolerance as we expect from using x=xO, 
meaning that xO is not a good approximation for the configuration at t + At. By solving now 
the 5 n q X 5 n q linear system JAx = — r using methods we will describe in next section, we 
correct xO towards a better solution. The while loop repeats the process, computes J and r 
again for the xO and x=xO+deltax and tests now the norm of the residual. If it passes, we 
assume x as being very close to ml, (/A and iCL that satisfy the residuals being equal to zero. 
The while loop exits and the NewtonRaphson function returns the output xO+deltax to the 
Update function. 

def Update(xO,hapO,k_lim,k_save,args) : 

#l.For a given time step k. 
for k in range (k_lim) : 

#2. Test whether current k is a multiple of k_step to save. 
if k°/„k_save==0 : 

"[Block of code] : save config" 
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#4- Set the first particular configurations as xO. 
x=array (xO) 

xO=Newton_Raphson (x , xO , hapO , hapO , args) 


which is redefined as the new xO. Notice above that all this procedure ocurred for a k=0 
iteration in the for loop, while now k=l starts with both xO and x as the previously obtained 
outputs of NewtonRaphson. 

The process itself is a sequence of k_lim NewtonRaphson algorithms each one composed 
by a sequence of solutions of linear systems of equations that ultimately yield a sequence of 
increments that correct each k iteration’s guess, x=xO, towards the pair (x,xO) that minimizes 
the residuals. 

By saving x at each k_save iteration, we have the magnetization and potential dynamics. 


5.6 A solver of linear systems of equations 

Each increment deltax is computed by the function rILU_AMG_GMRES which consists a Gen- 
eralized Minimum Residual Method (GMRES) we implemented from [25]. The algorithm is 
known to suffer from stagnation in the decrease of the residual and a standard solution is 
to accelarate the convergence by the use of Preconditioners [20]. Since we know the matrix 
structure we used for the non Poisson blocks the ILU preconditioner that was already avail- 
able in Python since we observed it was very fast. For the Poisson’s blocks of the Jacobian 
we implemented our version of the Algebraic Multigrid algorithm from [20, 23] and the coarse 
grid point selection algorithms from [26]. 
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Chapter 6 


Time Evolution of a Magnetic 
System 


The algorithms we developed so far for computing and assembling matrices, and our version of 
Newton- Raphson method that solves the JAx = — r, allow us to compute the time evolution of 
the system. We will compare our results for magnetization evolution and its Discrete Fourier 
Transform (DFT) with those for the standard problem given in [27], in order to evaluate the 
validity of our code and suggest improvements. 

6.1 Definition of the standard problem 

From [27], the system 0 consists of a flat box with shape 120 x 120 x 10 nm 3 . This is divided 
into cubic cells with the desired size, for example 5 x 5 x 10 nm 3 for Finite Differences where we 
associate a magnetization vector to each cell. For the Finite Element Method a 5 x 5 x 2.5 nm 3 
grainier mesh or a thinner mesh with 2.5 x 2.5 x 1 nm 3 is used that is then tessellated into 5 
or 6 tetrahedra [28] and to each vertex we associate a magnetization vector. 

The problem has now two stages: Relaxation stage and Dynamics stage. The first serves 
to define the initial configuration for the second. We start with a uniform magnetization 
m, = (0, 0, 1) T and an external applied magnetic field that is uniform across the system 
and independent of time given by h ap = (65.1, 46.5, 0) T kA/rn, which normalized by the 
x component yields h aPjn ps (1, 0.715, 0) T . Setting the damping constant to a high value, 
a c = 1, the evolution of each magnetization to align with the applied field is very fast and 
precession is supressed, as we want. We stop this stage once the equilibrium configuration is 
attained and save it for future use. It is suggested in [27] that 5 ns should suffice. 

The second stage then takes place, the goal is to record the dynamics of every magne- 
tization vector in P starting from the last configuration of the Relaxation stage. We reset 
the damping constant to the lower value of a c = 0.008 and redirect the applied field to 
h ap = (65.5, 45.9, 0) T kA/rn, also given by h app% « (1,0.7, 0) T , making a 35° with the x-axis. 
The update then runs for 20 ns with time step At and saving the magnetization configuration 
at every 5 ps. 

In both stages we use parameters for permalloy, the anisotropy of the system is not consid- 
ered, the saturation magnetization M s = 800 kA/rn, and exchange constant is 1.3 x 10 -11 J/m. 


67 



6.2 Data analysis 


The data obtained in the Dynamics stage is analysed using the Discrete Fourier Transform 
that we briefly summarize in Appendix K. 

We start by computing the average magnetization over the entire n q magnetization vectors 
in the sample as 

n q 

( m y) = —^2 m y,i ( 6 - 1 ) 

i= l 

and compare with the results from OOMMF micromagnetic simulator whose output is used 
as standard in [27]. 

Each magnetization vector in the z = 0 plane has its own evolution composed by N — 1 
steps, hence we have an N entry vector given as m y ^ = . . . ,m y) N- i)f where the 

set of all m.yfij was obtained from the Relaxation stage. It is on each one of these vectors 
that we apply the DFT to yield a vector of complex coefficients c,; as in (K.ll) where each 
entry is given by (K.12). There are N coefficients of these that will span frequencies from 
/o = 0 up to f N—i = (A" — 1)/(N At), where At is the time step of recorded data, 5 ns, hence 
we have fk/ GHz e [0, 200[. With the complex vectors, Cj, for all z = 0 vertices we choose one 
frequency, fk and average each amplitude as 

^ n xy 

(Ak) = V4, (6.2) 

n xv 

xy i = l 

where n xy is the number of magnetization vectors in that plane. The set of (Ak) is then plotted 
against fk- Again the results are compared with the same calculations for the magnetization 
evolution from OOMMF. 

We can go further now and analyse resonance peaks observed in this spectrum. By choos- 
ing their particular frequencies we plot the spatial distribution of the amplitudes and phases 
for each x.; in z = 0 slice of the mesh. 


6.3 Numerical Results of Relaxation Stage 

To mesh the system Q we used a BCC lattice. We chose a coarser mesh than suggested in [27] 
and used 5x5x5 nm 3 cells with a central vertex. Since the geometry is simple we developed 
our own mesh generator as well as the algorithms for computing the connectivity matrices 
where all distances were measured in units of Iex (c/ the fundamental scales in section 2.2). 

The initial configuration is given by the magnetization m.; = (0,0, l) 7 mentioned previ- 
ously as well as the the two potentials <fii and Ui at every vertex consistent with this magne- 
tizations field. Hence we have to solve first the coupled Poisson’s problem. 

After setting up the system for 5000 iterations with time step At = 1 ps measured in units 
of (7 e /^o Mg)^ 1 we simulated the 5 ns. The following evolution of magnetization y component 
was obtained 
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Figure 6.1: Evolution of the average magnetization components of (m) versus time. At 1.5ns we have 
approximately (.7886, .5913, —.0002) and at 5ns we have (.7884, .5892, —.0008) 


where we only show the first 1.5 ns since the magnetization is approximately constant for 
the remaining time. We note that the system seems to have reached an equilibrium state after 
0.5 ns. Further we note that the norm of magnetization evolves for all 5 ns as 




time (ns) 


Figure 6.2: Average norm and standard deviation versus time. 


meaning that the more iterations, the more we deviate from physically required norm 
conservation. A spatial plot shows that its the vertices at edges of system that most change 
their norm. The Dynamics stage will take longer than the first, since this issue will also emerge 
we decided to start it at a relaxed state but with few norms exceeding 1. Therefore instead 
of using the configuration at 5 ns, we decided to start the next stage with the configuration 
at 1.5 ns. The next figure shows a section with z = 0 of this data 
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Figure 6.3: Magnetization field for z — 0. The colorbar gives the z component. The grey magnetic field 
belongs to z — 0 plane and is the sum of the applied field with the much stronger stray field. The blue and 
red lines indicates the direction of the applied field for the relaxation and dynamics stages. 
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In the initial configuration all magnetization vectors point in the z direction and hence we 
had the four corners had the same local configuration. Upon introducing the applied field, 
the system evolved into a state where two corners on the left and right of the applied held 
have magnetization that follows the contour of the boundary and the bottom left and top 
right almost align with the applied held. 

The magnetization are mostly parallel hence the exchange held is small. Plot as a grey 
vector held is the sum of the applied held with the demagnetizing held, where the latter 
dominates. Observe the upper and right sides of the system we have a* > 0 while on the 
lower and left sides the charge density is negative. That is consistent with the direction of 
the demagnetizing held, it points in the direction opposite of the magnetization. Also notice 
that on the lower left and upper right sides the direction of h m is not what is expected. The 
demagnetizing held calculation on these two corners might not be accurate enough. 


6.4 Numerical Results of the Dynamical Stage 


We observed as our simulation proceeded that some of magnetization norms were not con- 
served. Additionally, since each NR iteration took from 10 s to 1 min, simulating 20 ns at 1 ps 
steps would require 20000 iterations. Hence we decided to use only data up to 5 ns. The m y 
results for the hrst 2 ns for our code and the OOMMF data from [27] are as follows 
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Figure 6.4: Evolution of magnetization y component for 2 ns for our simulation starting at 1.5 ns configuration 
of the Relaxation stage and OOMMF evolution starting at 5 ns. Initially: m y = .591; m® OMMF = .593 

The general shape of both waves is very similar. Our results start with magnetization, 
m y , 0.02 lower at t = 0, which can be caused by the progressive lowering of the magnetization 
norms as the number of iterations increased as shown in figure 6.2. The peaks are lower in 
our case and stretched being out of phase after 1.75 ns. The evolution of the y component 
of magnetization from 2 ns onward will keep oscillating with gradually smaller amplitudes 
but its average also becomes smaller, an aspect not observed in OOMMF whose oscillation is 
centered on m y = 0.586 for all 20 ns. 

The average amplitudes (6.2) of the DFT for each y component of each magnetization in 
the 2 = 0 plane show the following behaviour in a logarithmic base 10 scale 



Figure 6.5: Average amplitude for the y components for our simulation and OOMMF in logarithmic scale. 
The amplitude for k = 0 are 10~ 2 ' 7 for our results and 10~ 3 ' 3 for OOMMF. 

The OOMMF results use the full 20 ns giving denser data points than our 5 ns. It was 
expected from [27] that the bigger is the mesh cell size the higher are the amplitudes for 
the major peaks, which we do not observe in our results, our two major peaks are slightly 
lower than OOMMF. The absence of a denser number of data points smoothed out the main 


71 



peaks and those smaller ones for frequencies of 12.5 GHz onward. Additionally we observe 
the presence of very high amplitudes in lower frequencies in figure 6.5, that remains to be 
explained. 

We observe a close agreement for frequencies associated to the two main peaks, one at 
8.03 GHz and other at 10.91 GHz, which are shifted from the two peaks for OOMMF given 
respectively by 8.25 GHz and 11.25 GHz. The cause is associated with the cube dimensions 
used 5x5x5 nm 3 , and the Finite Element Method, as shown in [27] smaller mesh cube would 
lead to a better approximation to OOMMF results. The shift to the left of our spectrum 
shows that the major contributors for the evolution of the magnetization y component have 
smaller frequency, which is consistent with the magnetization evolution being stretched in 
figure 6.5. 

Now, choosing the resonance frequencies detected for the average magnetization we can 
analyse the spatial distribution of amplitudes and phases for the magnetization evolution in 
the z = 0 plane. 

From the frequency of the first resonance 8.03 GHz, we have the following amplitude 
distribution in logarithmic base 10 scale corrected by the minimum and the corresponding 
distribution of phases 



In the three components we observe for the amplitudes and phases the same symmetry 
already observed in figure 6.3. Additionally the phases are bigger at two corners coinciding 
with the smaller amplitudes. The same pattern is observed in the standard results in Appendix 
J. However figure 6.6 shows for the x and y components that amplitudes and phases at the 
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lower left and upper right corners are not as smooth as the standard results. This might be 
caused by the inaccurate stray field calculation as we observed in figure 6.3, according to [27] 
this field contribution is highly sensitive to the mesh structure. A thinner mesh structure 
could then improve the results or the code developed in section 4.5 must be improved at these 
boundary regions. 

The second peak at 10.91 GHz exhibits the following distributions of amplitudes and phases 



Figure 6.7: Spatial distribution of x , y and 2 amplitudes and phases for 10.91 GHz. 

Like the previous frequency we observe also good agreement between our results and the 
standard ones in figure J.2, with ours being less smooth. Additionally for the x component 
phases, we do not observe in figure 6.7 at the center of the system the 6 ~ n region we expect 
from the standard results. 

For the x component there is a yellow ” shell” enclosing a region with phases near — 1 while 
in z component we have — ir splitting the regions. This distribution shows domains of phases 
also observed in figure J.2. 

6.5 Conclusion 

Our results for the average magnetization evolution and DFT for the y component are in 
good agreement with the standard problem proposed in [27]. The (m y ) wave shows the 
same oscillations stretched which is consistent with the shift toward lower frequencies in the 
amplitude spectrum. 

The space resolved amplitude and phase are similar to those reported in [27] and shown in 
Appendix J. The non smooth amplitude and phases near the corners of the system coincide 
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with magnetization vectors with norm not conserved as the number of iterations increase as 
well as demagnetizing field having opposite direction to what is expected. Reviewing the code 
in section 4.5 at the system corners might solve the issue and/or using smaller cubic cells since 
demagnetizing fields will be more accurately computed. 

However the suggestion of reducing the mesh cell structure comes with the price of in- 
creasing the number of vertices in the mesh and then the size of the Jacobian matrix and 
residual vector, meaning larger computational times. 
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Chapter 7 


Conclusion and Outlook 


We started with micromagnetic theory that describes the dynamics of magnetization at the 
microscale. Two equations where introduced: Landau-Lisfshits equation and Landau-Lifshits- 
Gilbert equation. We rederived the first from the experimental observation that upon appli- 
cation of an external magnetic field magnetization precesses around this field and tends to 
align with it. From it we develop our intuition for the dynamics implied by each one of the 
four magnetic field contributions: the applied field; the anisotropy held that benefits certain 
directions within the system; the stray exchange held which tends to cancel hctitious magnetic 
charges and the exchange held which smooths out spatial variations in magnetization. 

Then we presented the LLG equation derived by Gilbert. Unlike LL, LLG implicitly 
determines the evolution of magnetization. We then established the goal of developing from 
scratch a Python code to integrate this equation. 

Defining the residuals for LLG and Poisson’s equation that determines the stray held, 
we used Finite Element Method to discretize space and Implicit Midpoint Rule to discretize 
time. Realizing that the stray held potential requires asymptotic boundary condition and to 
avoid the computational cost of having to handle the surrounding mesh we opted by a less 
expensive strategy. We introduced a new potential, solve its Poisson’s equation within the 
system and then use the Boundary Element Method equations to map this potential onto the 
values of the stray held potential at the boundary. As a consequence the previous asymp- 
totic boundary conditions were replaced by Dirichlet boundary conditions. This rendered 
integrating LLG equation together with two Poisson’s equation using only information about 
the system, avoiding the need for meshing the surrounding. To solve the equation we make a 
linear approximation as in Newton Raphson method. The linear system of equations obtained 
yields as a solution the update for a current conhguration of magnetization and potential in 
the system. 

Since this linear system has to be solved at every time step we developed the code to 
compute the matrix entries using as much as possible vectorized operations giving us fast 
calculations and assembly. 

To validate our code, we implemented the resolution of a standard problem, where we 
concluded that our results are close to those expected, and can be improved using a denser 
mesh. During these simulations we observed our solver of linear systems was the slowest step. 
We implemented our Python version of GMRES with two preconditioners, a Python already 
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available ILU preconditioner for non Poisson’s matrix blocks and our version of Algebraic 
Multigrid for the two Poisson’s blocks. Together they revealed to be at least five times faster 
than the Python sparse matrix solver, being absolute essential in this thesis. However much 
improvement needs to be done as the successive resolution of linear systems is what makes 
the code slow, not the actual assembly of matrix and residual vectors. Hence the main focus 
to improve has to be in using faster methods to solve the linear system of equations. 

Despite the high computational time involved, our code has as main feature its simplicity 
and clear relation with the numerical methods involved in micromagnetism: It serves the 
purpose of introducing the foundations of numerical micromagnetics, as a source of ideas and 
as a blueprint for more advanced codes. 
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Appendix A 


Data Structures for Magnetization, 
Potential and Gradients 


def fmalpha(m, alphae) : 

II II II 

m. shape=(3, nq) 

II II II 

ne = alphae . shape [1] 

malpha = zeros ( [3,4*ne] ) 

malpha[: , :ne] = m[:, alphae[0,:]] 

malpha [: ,ne:2*ne] = m[:, alphae [1,:]] 

malpha [: , 2*ne:3*ne] = m[:, alphae [2,:]] 

malpha [: , 3*ne:4*ne] = m[:, alphae [3,:]] 

return malpha 


Algorithm A.l: The magnetization at each (e, a) 

def f phialpha(phi , alphae) : 

II II II 

phi . shape=(nq, ) 

II II II 

ne = alphae . shape [1] 
phialpha = zeros ( [4*ne ,] ) 
phialpha[:ne] = phi [alphae [0] ] 

phialpha [ne : 2*ne] = phi [alphae [1] ] 

phialpha [2*ne:3*ne] = phi [alphae [2] ] 
phialpha [3*ne :4*ne] = phi [alphae [3] ] 
return phialpha 


Algorithm A. 2: The potential at each (e, a) 


def fmtot (malpha, ne) : 
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return (malpha [ : , : ne] +malpha [ : , ne : 2*ne] + 

malpha [ : , 2*ne : 3*ne] +malpha [ : , 3*ne : 4*ne] ) 


Algorithm A. 3: The total magnetization for each element e. 

def f gradient (ptot , alphae , signedvol , barcode) : 

ne = alphae . shape [1] #Number of elements=Number of cols. 

C = l/(6*signedvol) #1//X/. 

D_01 = ptot [: , alphae [0] ] -ptot [: , alphae [1] ] #Vector differences. 
D_02 = ptot [: , alphae [0] ] -ptot [: , alphae [2] ] 

D_03 = ptot [: , alphae [0] ] -ptot [: , alphae [3] ] 

D_12 = ptot [: , alphae [1] ] -ptot [: , alphae [2] ] 

D_13 = ptot [: , alphae [1] ] -ptot [:, alphae [3] ] 

G = zeros( [3,4*ne] ) #Initial G. 

G [ : , : ne] = array ( [ (D_12 [2 , : ] *D_13 [1 , : ] -D_12 [1 , : ] *D_13 [2 , : ] ) *C , 

(D_12 [0 , : ] *D_13 [2 , : ] -D_12 [2 , :]*D_13[0, :])*C, 

(D_12 [1 , : ] *D_13 [0 , : ] -D_12 [0 , :]*D_13[1, :])*C]) 

G [ : ,ne : 2*ne] = array ( [(D_02 [1 , :]*D_03[2, :]-D_02[2, :]*D_03[1, :])*C, 

(D_02 [2 , : ] *D_03 [0, :]-D_02[0, :]*D_03[2, :])*C, 
(D_02 [0 , : ] *D_03 [1 , : ] -D_02 [1 , : ] *D_03 [0 , : ] ) *C] ) 


G [ : , 2*ne : 3*ne] = array ( [(D_01 [2 , :]*D_03[1, :]-D_01[l, :]*D_03[2, :])*C, 

(D_01 [0 , : ] *D_03 [2 , : ] -D_01 [2 , : ] *D_03 [0 , : ] ) *C , 
(D_01 [1 , : ] *D_03 [0 , :]-D_01[0, :] *D_03[1 , :] )*C] ) 

G [ : , 3*ne:4*ne] = array ( [(D_01 [1 , :]*D_02[2, :]-D_01[2, :]*D_02[1, :])*C, 

(D_01 [2 , : ] *D_02 [0 , : ] -D_01 [0 , : ] *D_02 [2 , : ] ) *C , 
(D_01 [0 , : ] *D_02 [1 , : ] -D_01 [1 , : ] *D_02 [0 , : ] ) *C] ) 

with open_f ile ( ' Gradients°/oi ' °/ 0 (barcode) , 'w' ) as f: 
shape = G . shape 
atom = Float64Atom() 

Grads = f . create_carray(f . root Grads 1 , atom, shape) 

Grads [:] = G 

print ("File created & saved") 


Algorithm A. 4: The Gradients. Adapted from [24]. 
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Appendix B 


Poisson’s term 


Functions for the source term of Poisson’s residual for cf) and u 
def fdivm(malpha, Grads ,vol) : 

II II II 

IMPUT : To use IMR: malpha -> avgmalpha 
OUTPUT: An array with shape (ne,) with the 
divergences of magnetization at each element e. 

II II II 

ne = vol shape [0] 

#l.The gradients for each beta. 

GradsO = Grads [:,:ne] 

Gradsl = Grads [: ,ne : 2*ne] 

Grads2 = Grads [:, 2*ne : 3*ne] 

Grads3 = Grads [:, 3*ne :4*ne] 

#2. The magnetizations for each beta. 
malphaO = malpha ne] 
malphal = malpha [: ,ne : 2*ne] 
malpha2 = malpha [:, 2*ne : 3*ne] 
malpha3 = malpha [:, 3*ne : 4*ne] 

#3 .E-wise multiply , e-wise sum, sum column-wise to 
#give a row array with shape=(ne, ) . 

divm = sum(malphaO*GradsO+malphal*Gradsl+ 
malpha2*Grads2+malpha3*Grads3) 

return divm 


Algorithm B.l: Divergence of magnetization, (3.37). 
def fs_src_phi (malpha, vol, Grads, nq, barcode) : 

II II II 

The phi Poisson's residual source term contributes with 
this 5nqxl vector which will be modified in fs_phi_u function 
so that entries in blocks associated with a boundary vertex 


83 


are zero since we have the s=Gu-phi residual . 

OUTPUT: A csr matrix with shape (5*nq,l). 

II II II 

ne=vol . shape [0] 
s_src_phi = zeros ( [4*ne,] ) 

#1 .Divergence of magnetization for every e. 

divm = f divm(malpha, Grads, vol)*vol/4 
#2. Copy divm for every alpha 

s_src_phi [ :ne] = divm 

s_src_phi [ne : 2*ne] = divm 

s_src_phi [2*ne : 3*ne] = divm 
s_src_phi [3*ne :4*ne] = divm 

#3. Load indices and assemble the 5nqxl column vector. 

with open_f ile ( ' IgJg_s_src_phi“/„i ' % (barcode) , 'r' ) as f: 

Ig = f . root . IgJg_s_src_phi [ : , 0] 

Jg = f . root , IgJg_s_src_phi [ : , 1] 
s_src_phi = csr_matrix( (s_src_phi , (Ig, Jg) ) , shape=(5*nq, 1) ) 
return s_src_phi 


Algorithm B.2: The first term of Poisson’s residual, (3.39). 


def f s_src_phi_indices (alphae .barcode) : 

II II II 

Indices for the first term in Poisson's residual 
for phi. No indices for u are required since we 
will copy s_src_phi in fs_phi_u function for the 
u entries. 

OUTPUT: A file with Ig and Jg as the 1st 
and 2nd cols. 

II II II 

ne=alphae . shape [1] 

Ig = zeros ( [4*ne,] ,dtype =l int ' ) 

Jg = zeros ( [4*ne,] ,dtype =l int ' ) 

Ig[:ne] = alphae [0] *5+3 

Ig[ne:2*ne] = alphae [1] *5+3 
Ig [2*ne : 3*ne] = alphae [2] *5+3 
Ig [3*ne : 4*ne] = alphae [3] *5+3 

with open_f ile ( ' IgJg_s_src_phi°/„i ' "/(barcode), 'w') as f: 
shape=(Ig. shape [0] ,2) 
atom = Int32Atom() 

IgJg = f . create_carray(f root , 1 IgJg_s_src_phi ' .atom, shape) 

IgJg[: ,0] = Ig 
IgJg [ : > 1] = Jg 
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print ("s_src_phi indices saved.") 


Algorithm B.3: Indices to map each s_divm entry into the residual vector. 


Functions for the stiffness term of Poisson’s residual for </> and u 


def JPoisson(vol ,ne , Grads, barcode) : 

II II II 

Stiffness matrix entries for Poisson's residual and Jacobian. 
For the first IMR will not require Kg*0.5 as we multiply 
Kg matrix by the average configuration. If we use Kg 
in a Jacobian we have to multiply Kg*0.5. 

OUTPUT: A file with the array (16*ne,) with all entries. 

II II II 


G_0 = Grads [:,:ne] 

G_1 = Grads [: ,ne : 2*ne] 

G_2 = Grads [: ,2*ne : 3*ne] 

G_3 = Grads [: ,3*ne : 4*ne] 

Kg = zeros( [ne*16,] ) 

#(alpha, beta) 

Kg[0:ne] 

= sum(G_0*G_0) *vol 

#(0,0) 

Kg [ne : 2*ne] 

= sum(G_0*G_l) *vol 

#(0,1) 

Kg [2*ne : 3*ne] 

= sum(G_0*G_2) *vol 

#(0,2) 

Kg [3*ne : 4*ne] 

= sum(G_0*G_3) *vol 

#(0,3) 

Kg [4*ne : 5*ne] 

= Kg [ne : 2*ne] 

#(1,0) 

Kg [5*ne : 6*ne] 

= sum(G_l*G_l) *vol 

#(1,1) 

Kg [6*ne : 7*ne] 

= sum(G_l*G_2) *vol 

#(1,2) 

Kg [7*ne : 8*ne] 

= sum(G_l*G_3) *vol 

#(1,3) 

Kg [8*ne : 9*ne] 

= Kg [2*ne : 3*ne] 

#(2,0) 

Kg [9*ne : 10*ne] 

= Kg [6*ne : 7*ne] 

#(2,1) 

Kg[10*ne : ll*ne] 

= sum(G_2*G_2) *vol 

#(2,2) 

Kg [ll*ne : 12*ne] 

= sum(G_2*G_3) *vol 

#(2,3) 

Kg [12*ne : 13*ne] 

= Kg[3*ne:4*ne] 

#(3,0) 

Kg[13*ne : 14*ne] 

= Kg [7*ne : 8*ne] 

#(3,1) 

Kg [14*ne : 15*ne] 

= Kg [ll*ne : 12*ne] 

#(3,2) 

Kg [15*ne : 16*ne] 

= sum(G_3*G_3) *vol 

#(3,3) 


with open_f ile ( ' PoissonJacobian°/„i ' % (barcode), 'w') as f: #Save 
shape=Kg . shape 
atom = Float64Atom() 

KgPoisson = f . create_carray (f root, 'Kg' , atom, shape) 
KgPoisson[:] = Kg 
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print("File created & saved") 


return Kg 


Algorithm B.4: Stiffness matrix entries, (3.44). Adapter from [24]. 


Function for the boundary term of Poisson’s residual for the potential u 
def fs_u_BC(malphas, alphas, normal_vecs, areas, nq) : 

II II II 

Poisson's residual surface integral for the potential u. 

Similar algorithm to BC_Q_u function. 

INPUT: 

malphas is analogous to malpha but comes from using instead alphas, 
malpha. shape=(3 , 3*ns) 

To use IMR: malphas -> avgmalphas 

OUTPUT: A csr array with shape (5*nq,l). Out of the 
three terms for Poisson's residual , the surface integral 
has a minus sign which is not included in this function. 

Hence we sum the 1st and 2nd terms and subtract this one: 
sphiu=s_divm+Kg . dot (x) -s_u_BC in fs_phi_u function. 

II II II 

ns=alphas . shape [1] 

BC=zeros ( [3*3*ns , ] ) 

Ig=zeros ( [3*3*ns , ] ) 

Jg=zeros ( [3*3*ns , ] ) 

C=areas/ 12 
i=0 

kd = lambda a,b: 1 if a==b else 0 #Kronecker Delta fc. 
for b in range (3): 

mbeta = malphas [: ,b*ns : (b+1) *ns] 
mbetan = sum(m_b*normal_vecs) *C 
for a in range(3): 

mbetandelta=mbetan* (l+kd(a,b) ) 

BC [i*ns : (i+1) *ns] =mbetandelta 
Ig [i*ns : (i+1) *ns] =alphas [a] *5+4 
i+=l 

return csr_matrix( (BC, (Ig, Jg) ) , shape=(5*nq, 1) ) 


Algorithm B.5: Third contribution for Poisson’s residual, (3.52). 
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Residual for Poisson’s equation for cj) and u potential plus the s^ = Guj nt — 4>i nt 
boundary term 

def f s_phi_u (KgphiuG , x , malpha , malphas , vol , Grads , Igs , 
nq, alphas ,normal_vecs , areas .barcode) : 

// II II 

Poisson's residual for phi and u plus s=Gu-phi . 

INPUT: Kg_phi_ [-1] = Stiffness matrix for phi with 
boundary entries replaced by -1 . 

KgphiuG -> Kg_phi_ [-IJ + G + Kg_u. 
x is the current configuration , shape=(5*nq, ) . 

To use IMR: x, malpha and malphas -> averages . 
malphas is analogous to malpha but uses alphas in its 
construction. 

OUTPUT: A (5*nq,l) csr vector with all Poisson's residuals . 

II II II 

#l.The source term. 

s_src_phi_u=f s_src_phi (malpha, vol , Grads ,nq, barcode) 

s_src_phi_u [4 :: 5] =s_src_phi_u [3 :: 5] #copy src term of phi to u entries. 
s_src_phi_u[Igs*5+3] =0 #zero the phi boundary entries. 

#2. The boundary u term. 

s_u_BC=f s_u_BC (malphas , alphas , normal_vecs , areas , nq) 

#3 . Compute the residual (KgphiuG does not have 0.5, its on avg x) 

sphiu=s_src_phi_u+KgphiuG dot (x) -s_u_BC 
return sphiu 


Algorithm B.6: Computes the residual for both potentials given by (4.56) for vertices within the 
system, at boundary by (4.60) and for system plus boundary by (4.64). 
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Appendix C 


Exchange Term 


def JEx(malpha,mtot .barcode ,nq,ne) : 

II II II 

Produces BO and B flattened. The first 

is used for rexchange, the second in the Jacobian. 

INPUT: malpha -> avgmalpha; mtot -> avgmtot. 

II II II 

ii=outer ( [0, 1 ,2 ,3] , [1 , 1 , 1 , 1] ) . f lattenO #Sequence of alphas (ii) 
j j=outer ( [1 , 1 , 1 , 1] , [0 , 1 , 2,3] ) . f lattenO #and betas (jj). 

#l.#Load K matrix. 

with open_f ile ( ' PoissonJacobian°/„i 1 % (barcode), 'r') as f: 

Kg=f , root . Kg [ : 16*ne] 

#2. Compute first term. 

BO = zeros( [3,4*ne] ) 
for c in range(16): 

a=ii[c] ttalpha. 

b=jj [c] #beta. 

BO [ : , a*ne : (a+1) *ne] +=malpha [ : ,b*ne : (b+1) *ne] *Kg [c*ne : (c+1) *ne] 
#3 . Compute second term and add to the first. 

Bvec = zeros ( [3 , 16*ne] ) 
for c in range(16): 

a=ii[c] #alpha. 

Bvec [ : , c*ne : (c+1) *ne] = (B0 [ : , a*ne : (a+1) *ne] - 

Kg [c*ne : (c+1) *ne] *mtot) /4 

ttf .Flatten twice. 

B=zeros( [16*ne*3*2,] ) 

#Lower triangular . Flatten and apply (-1) to y components 
B[:16*ne] = Bvec[0,:] #x 

B [16+ne : 2*16*ne] = -Bvec[l,:] tty 

B [2*16*ne : 3*16*ne] = Bvec [2,:] ttz 
ttUpper triangular has oposite signal. 

B[3*16*ne:] = -B[:3*16*ne] 
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#5. Save. 

with open_f ile ( ' IgJgJEx°/„i 1 % (barcode), 'r') as f: 
Ig=f , root IgJg JEx [ : , 0] 

Jg=f.root IgJgJEx[:,l] 

B=csr_matrix( (B , (Ig, Jg) ) , shape=(5*nq, 5*nq) ) 
return [B0,B] 


Algorithm C.l: Exchange contribution for the Jacobian, (3.93). 

def JEx_indices (alphae .barcode) : 
ne=alphae , shape [1] 

#l.For x row of Bvec build the sequence of 
#(alpha (i) , gamma (j)) indices 

ii=outer ( [0, 1,2,3] , [1,1, 1,1]) . f lattenO 
j j=outer ( [1,1, 1,1] , [0, 1,2,3]) . f lattenO 
#2. For a row of Bvec build the (16) block indices 
Ig = zeros ( [16*ne,] ,dtype =l int ' ) 

Jg = zeros ( [16*ne,] ,dtype =l int ' ) 
for c in range(16): 

Ig [c*ne : (c+1) *ne] =alphae [ii [c] , : ] 

Jg [c*ne : (c+1) *ne] =alphae [j j [c] , : ] 

#3. Build the arrays to receive global indices of every entry in B. 
#The structure of B determines the shape of IgExch and JgExch: 

#3 -> x,y,z, 16*ne -> row lenght of Bvec, 

#2 -> the lower triangular part of fxf blocks. 

IgEx=zeros( [3*16*ne*2,] ,dtype= 1 int 1 ) 

JgEx=zeros( [3*16*ne*2,] ,dtype= 1 int 1 ) 

#f.To the 16 block locations in Jacobian, 

#add the local x,y,z increments. 

Igxyzf ix=array ([2,2,1]) 

Jgxyzf ix=array ([1,0,0]) 
for d in range (3): 

#lower trianglar. 

IgEx [d*ne*16 : (d+1) *ne*16] =Ig*5+Igxyzf ix [d] 

JgEx [d*ne*16 : (d+1) *ne*16] = Jg*5+ Jgxyzf ix [d] 

#upper triangular. 

IgEx [(3+d) *ne*16 : (3+d+l) *ne*16] =Ig*5+ Jgxyzf ix [d] 

JgEx [(3+d) *ne*16 : (3+d+l) *ne*16] =Jg*5+Igxyzf ix [d] 

#f.Save. 

with open_f ile ( ' IgJgJEx°/„i 1 °/ 0 (barcode) , 'w' ) as f: 
shape= (IgEx . shape [0] ,2) 

atom = Int32Atom() #lst col: Ig; 2nd col: Jg 

IgJgJEx = f . create_carray(f . root , 1 IgJgJEx ', atom, shape) 

IgJgJEx [ : , 0] = IgEx 
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IgJgJEx [ : , 1] = JgEx 
return IgEx , JgEx 


Algorithm C.2: Indices for C.l. 


def rEx(B0 ,mtot ,nq,ne , barcode) : 

// II II 

INPUT: mtot -> avgmtot; BO is computed by JEx function. 

OUTPUT: a csr matrix with shape=(5*nq, 1) of exch residuals . 

II II II 

#1. Cross is computed for each vertex alpha 
rvec=zeros( [3,4*ne] ) 
for a in range (4): 

x = mtot [ [1 ,2 , 0] , : ] *B0 [ [2 ,0 , 1] , a*ne : (a+1) *ne] #Permute rows in mtot 
y = mtot [ [2 ,0 , 1] , : ] *B0 [ [1 ,2 , 0] , a*ne : (a+1) *ne] #and BO and e-wise mul:. 
rvec [ : , a*ne : (a+1) *ne] = -(x-y)/4 #-(cross product )/f.. 

#2. Flatning 
r_EX=zeros ( [3*4*ne , ] ) 
r_EX[:4*ne] = rvec [0] 

r_EX [4*ne : 8*ne] = rvec[l] 
r_EX[8*ne:] = rvec [2] 

#3. Save 

with open_f ile ( ' IgJgExres°/ 0 i ' °/ 0 (barcode) , 'r') as f: 

Ig_EX = f root . IgJgExres [:, 0] 

Jg_EX = f root . IgJgExres [:, 1] 
return csr_matrix( (r_EX, (Ig_EX, Jg_EX) ) , shape= (5*nq, 1) ) 


Algorithm C.3: Exchange contribution for the LLG residual, 
def rEx_indices (alphae , barcode) : 

II II II 

Builds and saves all indices for exchange residual . 

The file is then loaded by rEx function. 

The indices are also the ones required by the rD term. 

II II II 

#1. Define indices for the x component. 
ne=alphae . shape [1] 

Ig=zeros ( [4*ne , ] , dtype= 1 int ' ) 
for a in range (4): 

Ig [a*ne : (a+1) *ne] =alphae [a] 

#2. Build indices for each one of the three rows of 
#rvec in rEx function. 

Igxyzf ix=array ( [0 , 1 , 2] ) #x, y ,z componets 
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IgEx=zeros ( [3*4*ne , ] ) 
for d in range (3): 

IgEx [d*4*ne : (d+1) *4*ne] =Ig*5+Igxyzf ix [d] 

#3. Save the indices. 

with open_f ile ( ' IgJgExres°/oi ' °/ 0 (barcode) , 'w') as f: 
shape= (IgEx . shape [0] ,2) 

atom = Int32Atom() #lst col: IgEx; 2nd col: JgEx. 

IgJgExres = f . create_carray (f root , 1 IgJgExres atom, shape) 
IgJgExres [ : , 0] =IgEx 

IgJgExres [ : , 1] =zeros ( [3*4*ne , ] , dtype= 1 int ' ) 
print ( "rEx_indices saved, (also used for rD)") 


Algorithm C.4: Indices for C.3. 
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Appendix D 


Applied field, Stray Field and 
Damping 


The strategy to compute this Jacobian is to put the three term into arrays with shape 3 x 4n e , 
four blocks associated with a and n e vector columns each. We start by element-wise multiply 
each one of the four blocks of phialpha with the corresponding block from Grads and add 
them resulting in a 3 x n e array, which is then added to each block of happalpha and malpha 
yielding a 3 x 4n e array Dvec with four a blocks. It is then flattened twice into a single row 
with shape 1 X (2*3*4* n e ), two halves: one for lower triangular and other with opposite 
sign for upper triangular. The indices IgD and JgD are computed starting with the Ig indices 
for a single Dvec row. Since the blocks all belong to the main J diagonal all we need is Ig. 
We use it repeatedly to build indices for the six blocks of D by multiplying all of them with a 
factor of 5 and then adding the right increments of +0, +1 and +2 depending on whether its 
a block associated with x, y or z components on the lower or upper triangular parts. 

The residual calculation uses similar strategies as JD, the new aspect is the calculation of 
the cross product already explained in Chapter 5. Observe the indices for rD coincide with 
those for rEx. 

def JD (happalpha , phialpha , delt at , Grads , malphaO , vol , nq , ne , alphac , barcode) : 

II II II 

Builds all entries for the D term. 

INPUT: happalpha and phialpha -> averages 

except malphaO; Grads is the gradient matrix with shape=(3,4*ne) ; 

vol is the positive volumes of all tetrahedra; 

alphac is the damping constant ; barcode identifies the file 

with the indices IgD and JgD to be loaded. 

OUTPUT: A 5*ng by 5*nq csr matrix. 

II II II 

#1. Define constants. 

vol=vol/4 

cte=alphac/deltat 

#2. Gradient of phi for each element e. 

gradphialpha=zeros ( [3 , ne] ) 
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for a in range (4): 

gradphialpha+=Grads [ : ,a*ne: (a+l)*ne] *phialpha[a*ne : (a+l)*ne] 

#3.D vectors to be mapped, in skew matrix. 

Dvec=zeros( [3,4*ne] ) 
for a in range (4): 

Dvec [ : , a*ne : (a+1) *ne] =(-happalpha [ : , a*ne : (a+1) *ne] *0 . 5 

+gradphialpha*0 . 5 

-cte*malphaO [ : , a*ne : (a+1) *ne] ) *vol 
ttf. Flatten Dvec to get the lower triangular part and 
#copy* (-1) for the upper triangular. 

D = zeros ( [2*3*4*ne ,] ) 

#Lower triangular . 

D[:4*ne] = Dvec[0] #x 

D[4*ne:8*ne] =-Dvec[l] tty 
D [8*ne : 12*ne] = Dvec [2] ttz 
ttUpper triangular . 

D [12*ne : ] =-D [ : 12*ne] 

#5. Load indices and build the matrix in csr format. 

with open_f ile ( ' IgJgJD°/ 0 i 1 "/(barcode) , 'r 1 ) as f: 

IgD=f . root , Ig Jg JD [ : , 0] 

JgD=f . root . Ig Jg JD [ : , 1] 

return csr_matrix( (D , (IgD, JgD)) , shape=(5*nq,5*nq) ) 


Algorithm D.l: Applied field, stray field and damping Jacobian. 

def JD_indices (alphae ,ne .barcode) : 

II II II 

All the indices ig,jg for the D term entries 

computed with JD function and saved at 'IgJgJD/i' ‘/(barcode) . 

II II II 

ttl. Define indices for the x component. 

Ig=zeros ( [4*ne , ] ) 
for a in range (4): 

Ig [a*ne : (a+1) *ne] =alphae [a] 

#2. Build indices for other components and final 5*nqx5*nq Jacobian. 

Igxyzf ix=array ([2,2,1]) 

Jgxyzf ix=array ( [1,0,0]) 

IgD=zeros ( [2*3*4*ne , ] ) 

JgD=zeros ( [2*3*4*ne , ] ) 
for d in range (3): 

ttLower triangular 

IgD [d*4*ne : (d+1) *4*ne] =Ig*5+Igxyzf ix [d] 

JgD [d*4*ne : (d+1) *4*ne] =Ig*5+Jgxyzf ix [d] 

ttUpper triangular 


94 


IgD [(3+d) *4*ne : (3+d+l) *4*ne] =Ig*5+Jgxyzf ix [d] 

JgD [(3+d) *4*ne : (3+d+l) *4*ne] =Ig*5+Igxyzf ix [d] 

#3. Save the indices. 

with open_f ile ( ' IgJgJD°/ 0 i 1 % (barcode) , 'w 1 ) as f: 
shape= (IgD shape [0] ,2) 
atom=Int32Atom() 

IgJgJD = f create_carray (f , root , 1 IgJgJD atom, shape) 
IgJgJD [ : > 0] = IgD 
IgJgJD [: ,1] = JgD 


Algorithm D.2: Indices for D.l. 

def rD (happalpha, Grads ,phialpha,malpha,dmalpha, 
alphae , vol , nq , ne , deltat , alphac , barcode) : 

II II II 

Computes the residual vector for the applied field+stray field+ 
damping terms of LLG. 

INPUT: happalpha -> av ghapp alpha ; malpha -> avgmalpha; 
phialpha -> avgphialpha; dmalpha = malpha-malphaO (not avgmalpha) 
OUPUT:A 5nqxl residual vector is csr format. 

II II II 

#1 . Constants . 

rvec=zeros( [3,4*ne] ) 
vol=vol/4 
cte=alphac/ deltat 

#2. Gradient of phi for each element e. 

gradphialpha=zeros ( [3 , ne] ) 
for a in range (4): 

gradphialpha+=Grads [ : ,a*ne: (a+l)*ne] *phialpha[a*ne : (a+l)*ne] 

#3. Vectors to be mapped into the residual vector. 

Dvec=zeros( [3,4*ne] ) 
for a in range (4): 

Dvec [ : , a*ne : (a+1) *ne] = (+happalpha [ : , a*ne : (a+1) *ne] 

-gradphialpha 

-cte*dmalpha [ : , a*ne : (a+1) *ne] ) *vol 

#f. Cross product. 
for a in range (4): 

malphal20=malpha[ [1,2,0] ,a*ne: (a+l)*ne] 
malpha201=malpha[ [2 ,0, 1] ,a*ne : (a+1) *ne] 
x=malphal20*Dvec [[2,0,1] ,a*ne : (a+1) *ne] 
y=malpha201*Dvec [[1,2,0] ,a*ne : (a+1) *ne] 
rvec [ : , a*ne : (a+1) *ne] =x-y 
#5 . Flatning . 
r_D=zeros ( [3*4*ne , ] ) 
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r_D[:4*ne] = rvec [0] 

r_D [4*ne : 8*ne] = rvec[l] 
r_D[8*ne:] = rvec [2] 

#6. Build the residual vector using IgrD and JgrD, 
#that are the same as IgExchres and JgExchres . 

with open_f ile ( ' IgJgExres°/ 0 i ' % (barcode) , 1 r ' ) as f: 
IgrD=f , root , Ig JgExres [ : , 0] 

JgrD=f . root . IgJgExres [ : , 1] 
return csr_matrix( (r_D, (IgrD, JgrD) ) , shape=(5*nq, 1) ) 


Algorithm D.3: Applied field, stray field and damping contributions for the LLG residual. 
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Appendix E 


P and Q terms 


Both terms are computed by cycling through all combinations of (a, /3) and evaluating the 
3x1 expressions for (3.77) or the first term in (3.56) giving us Pvec or Qvec arrays with shape 
3 x 16n e . These arrays are then flattened giving us P and Q 1 x 3 * 16n e arrays. To map these 
entries into the correct locations within the 5x5 block (5.3) for all blocks in J we build the 
indices as follows. Observing P and Q arrays has a (a, j3) 16 entry sequence of indices three 
consecutive times, each one requiring a local fix corresponding to x, y and z components of 
the vectors. Since both have the same flattened structure, observing in (5.3) the locations of 

and are the transpose of one to the other, we conclude that the local transformations 
given by IgPxyzfix[d] and JgPxyzfix[d] have to switch roles when obtaining the indices 
of Q from Ig and Jg. 

Since LLG does not depend on the potential u we do not have a P term on the fifth column 
of a 5 x 5 block. However, the Poisson’s equation for the potential u has a Q term that is 
equal to the one for but is located one row below, hence all we have to do is use JgP+1 and 
IgP. 

Observe also from the Section 4.6 that at every row block 5 i x 5 n q that is associated to a 
vertex i e dQ we do not have a Poisson’s residual but (4.60), hence the entries for these rows 
of Q for the potential f> have to set equal to zero. 

Finally IMR is not implemented in either of JP_phi or JQ_phi, the malpha input has to 
be equal to (malpha+malphaO) *0 . 5. 

ds 

For the boundary vertices the entries for have an additional term given by the second 
term in (3.56) which is computed with BC_Q_u. 

def JP_phi (malpha, Grads ,nq,vol .barcode) : 

II II II 

INPUT: malpha -> avgmalpha 

OUTPUT: The P contribution for the Jacobian in csr format. 

To include IMR use avgmalpha and multiply by an extra 1/2 

from the potential . 

The 1/8 factor in vol comes from: 1/2 from avgphi and 

1/f from basis integral . 

II II II 

#1 . Constants . 
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ii = outer ( [0, 1 ,2 ,3] , [1 , 1 , 1 , 1] ) . f lattenO 
jj = outer ( [1 , 1 , 1 , 1] , [0 , 1 , 2,3] ). f lattenO 
vol= vol/8 
ne = vol shape [0] 

#2. Compute the cross procucts between m_a and grad(N_b) . 
Pvec=zeros ( [3 , 16*ne] ) 
for c in range(16): 

malphal20 = malpha[ [1 , 2, 0] , ii [c] *ne : (ii [c] +1) *ne] 
malpha201 = malpha[ [2, 0, 1] , ii [c] *ne : (ii [c] +1) *ne] 
a = malphal20*Grads [ [2,0 , 1] , j j [c] *ne : ( j j [c] +1) *ne] 
b = malpha201*Grads [ [1 ,2 ,0] , j j [c] *ne : ( j j [c] +1) *ne] 

Pvec [ : , c*ne : (c+1) *ne] =- (a-b) *vol 
#3 . Flatten. 

P=zeros( [3*16*ne,] ) 

P[:16*ne] = Pvec [0] #x 

P [16*ne : 32*ne] = Pvec[l] tty 

P [32*ne : 48*ne] = Pvec [2] ttz 

.Assemble the 5*nq by 5*nq P matrix. 
with open_f ile ( ' IgJgJPQindices_phi%i 1 °/ 0 (barcode) , 'r') as f: 
IgP = f . root IgJgJPQ_phi [: ,0] 

JgP = f.root IgJgJPQ_phi [ : , 1] 
return csr_matrix( (P , (IgP, JgP)) , shape=(5*nq,5*nq) ) 


Algorithm E.l: P contribution for the Jacobian, (3.77). 
def JQ_phi (Grads , vol) : 

II II II 

Computes all entries for the Q term without taking 
into account the surface integral contribuition. 

Notice the output has no minus sign. 

Unlike JP, this term is independent of the configuration. 

To include IMR multiply the output by 1/2. 

OUTPUT: A (3*16*ne, ) array with all the entries. 

II II II 

#1 . Constants . 

vol=vol/4 
ne=vol shape [0] 

j j=outer ( [1,1, 1,1] , [0, 1,2,3]) . f lattenO 

#2. For each alpha cycle through all beta's and 

#arrange grad(N_b) in that sequence. 

Qvec=zeros ( [3 , 16*ne] ) 
for c in range(16): 

Qvec [ : , c*ne : (c+1) *ne] =Grads [ : , j j [c] *ne : (j j [c] +1) *ne] *vol 

#3 .Flatten. 
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Q=zeros( [3*16*ne,] ) 
Q[:16*ne] = Qvec [0] #x 

Q [16*ne : 32*ne] = Qvec[l] #y 
Q [32*ne : 48*ne] = Qvec [2] #z 
return Q 


Algorithm E.2: First part of Q term, (3.56), for both potential <f> and u 

def JPQ_indices_phi (alphae , barcode) : 

II II II 

The indices for P are IgP, JgP. The 

indices for Q_phi are IgQ, JgQ. For first term 

of Q_u use IgQ+1, JgQ. (there no P column for u) 

0UTPUT:A file with an array where the 1st col 
is IgP and the 2nd is JgP, the 3rd is IgQ and 
the ^th is JgQ. All with shape=(3*16*ne, ) . 

II II II 

#1 . Compute the 16 sequence of alphas and betas 

ne=alphae . shape [1] 

ii = outer ( [0, 1 ,2 ,3] , [1 , 1 , 1 , 1] ) . f lattenO 
jj = outer ( [1 , 1 , 1 , 1] , [0 , 1 , 2,3] ). f lattenO 
#2. The global vertices indices for one row of Pvec. 

Ig = zeros ( [16*ne,] ,dtype =l int ' ) 

Jg = zeros ( [16*ne,] ,dtype=' int 0 
for c in range(16): 

Ig [c*ne : (c+1) *ne] = alphae [ii [c] , :] 

Jg [c*ne : (c+1) *ne] = alphae [jj [c] , :] 

#3. From Ig and Jg fix them to locate the 
#vectorial contribution of P in the 5x5 blocks. 

IgP=zeros( [3*16*ne,] ,dtype=' int 1 ) 

JgP=zeros( [3*16*ne ,] ,dtype= ' int 1 ) 

IgQ=zeros( [3*16*ne,] ,dtype=' int 1 ) 

JgQ=zeros( [3*16*ne ,] ,dtype= ' int 1 ) 

IgPxyzf ix=array ( [0, 1 , 2] ) #x, y, z 

JgPxyzf ix=array ( [3,3,3] ) #fouth column of 5x5 block. 
for d in range (3): 

IgP [d*16*ne : (d+1) *16*ne] =Ig*5+IgPxyzf ix [d] #Local fix 
JgP [d*16*ne : (d+1) *16*ne] =Jg*5+JgPxyzf ix [d] #for P indices. 

IgQ [d*16*ne : (d+1) *16*ne] =Ig*5+JgPxyzf ix [d] #For Q, switch 
JgQ [d*16*ne : (d+1) *16*ne] =Jg*5+IgPxyzf ix [d] #local increments. 
Save indices. 

with open_f ile ( ' IgJgJPQ_phi°/oi 1 % (barcode) ,' w ' ) as f: 
shape=(IgP shape[0],4) 

atom = Int32Atom() #lst col: IgP; 2nd col: JgP 
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IgJgJPQ = f , create_carray(f . root ,' IgJgJPQ_phi ', atom, shape) 

IgJgJPQ L : > 0] = IgP 
IgJgJPQ [:,1] = JgP 
IgJgJPQ L : > 2] = IgQ 
IgJgJPQ [: >3] = JgQ 
print ( "JPQ_indices saved.") 


Algorithm E.3: Indices for both P and Q terms for </> potential, for the u potential there is only the Q 
term. 


def BC_Q_u(normal_vecs , alphas , areas ,nq) : 

II II II 

For every (a,b) comb there is a 3xns array of normal 
vector that have to be mapped. Chosing (a,b) we 
flatten normal_vecs by picking one of its rows 
and build the corresponding indices. 

To include IMR multiply the output by 1/2. 

Similar algorithm to fs_u_BC. 

OUTPUT: A 5*nqx5*nq csr array. 

II II II 

ns=alphas , shape [1] 

Q_BC=zeros ( [ns*3* (3*3) , ] ) 

Ig=zeros ( [ns*3* (3*3) , ] ) 

Jg=zeros ( [ns*3* (3*3) , ] ) 

C=(areas/12) 

i=0 

kd = lambda a,b: 1 if a==b else 0 #Kronecker Delta fc. 
for a in range (3): 

for b in range (3): 

for x in range (3): 

Q_BC [i*ns : (i+1) *ns] =-normal_vecs [x] * (l+kd(a,b) ) *C 
Ig [i*ns : (i+1) *ns] =alphas [a] *5+4 
Jg [i*ns : (i+1) *ns] =alphas [b] *5+x 
i+=l 

return csr_matrix( (Q_BC, (Ig, Jg)) , shape=(5*nq,5*nq) ) 


Algorithm E.4: The second term in (3.56) for the u potential. 
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Appendix F 


Time derivative term 


The time derivative Jacobian is build starting with a 3 x 4n e array, Avec, where each a block 
has n e equal column vectors given by 4A ) (1. 1, 1) T . Then we flatten for each component giving 
us A and build the corresponding indices that map these entries to the main diagonal of the 
(i,j) block (3.75). 

The residual only requires multiplying each a block of dmalpha=malpha-malphaO by the 
n e volumes. The final structure is the same as the exchange residual contribution. 

def fA(alphae,vol,deltat) : 
vol=vol/ (deltat*4) 

Avec=ones ( [3 , 4*ne] ) 

#E-wise multiply each 3xne block of Avec with volumes. 
for a in range (4): 

Avec [ : , a*ne : (a+1) *ne] *=vol 
#Flatten . 

A=zeros ( [3*4*ne , ] ) 

A[:4*ne] = Avec[0] 

A [4*ne : 8*ne] = Avec[l] 

A [8*ne : 12*ne] = Avec [2] 
return A 


Algorithm F.l: Entries given by (3.64). 

def JA_indices (alphae , barcode) : 
ne=alphae . shape [1] 

IgA=zeros ( [3*4*ne , ] , dtype= 1 int ' ) 

for a in range(4): #diagonal entry: 

IgA [a*ne : (a+1) *ne] = alphae [a] *5 #lst 

IgA[4*ne+a*ne:4*ne+(a+l)*ne] = alphae [a] *5+1 #2nd 
IgA [8*ne+a*ne : 8*ne+(a+l) *ne] = alphae [a] *5+2 #3rd 
with open_f ile ( ' IgJgA°/ 0 i 1 % (barcode) , ' w 1 ) as f: 
shape=(IgA shape[0],l) 
atom=Int32Atom() 
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IgJgJA=f create_carray (f root, 'IgJgJA' , atom, shape) 

IgJgJA [ : , 0] =IgA 

print ( " JA_indices saved . " ) 


Algorithm F.2: Indices for F.l. 


def rA(dmalpha,deltat,nq,ne,vol, barcode) : 
vol = vol/ (4*deltat) 

#Multiply each alpha block by the volumes/ (/*deltat) 
for a in range (4): 

dmalpha [ : ,a*ne : (a+1) *ne] *=vol 
#Flatten. 

r_A = zeros ( [4*ne*3,] ) 
r_A[:4*ne] = dmalpha [0 , : ] 

r_A [4*ne : 8*ne] = dmalpha [1,:] 
r_A[8*ne:] = dmalpha [2,:] 

with open_f ile ( ' IgJgExres°/ 0 i ' °/ 0 (barcode) , 1 r ' ) as f: 
Ig_A = f root IgJgExres [: ,0] 

Jg_A = f . root IgJgExres [:, 1] 
return csr_matrix( (r_A, (Ig_A, Jg_A) ) , shape=(5*nq, 1) ) 


Algorithm F.3: Residual for the time derivative term, (3.63) 


102 


Appendix G 


Assembly J and r 


def J_fix_assembly (barcode) : 

#O.Load necessary data. 

#Connectivity matrix. 

with open_f ile ( ' alphae°/ 0 i 1 % (barcode) , 'r 1 ) as f: 

alphae = f root . alphae [ : , : ] 
ne = alphae . shape [1] 

#Connectivity array for boundary vertices . 

#Vetices at boundary are Igs. 

with open_f ile ( ' alphas°/„i ' % (barcode) , 'r 1 ) as f: 
alphas = f root . alphas [ : ] 

Igs = f root Igs[:] 
ns = alphas . shape [1] 

#Normal vectors. 

with open_f ile( 'normalvecs7 0 i ' % (barcode) , 1 r ' ) as f: 
normal_vecs = f root normalvecs [ : ] 

#Vertices positions . 

with open_f ile ( ' allvertexlocations°/oi ' 7« (barcode) r ' ) as f: 

ptot = array (f . root .ptot [:] ) 
ptot = transpose (ptot) 
nq = ptot . shape [1] 

#Poisson Jacobian for phi. For u potential add +1 . 

with open_f ile ( ' IgJgJPoisson_phi%i ' °/„ (barcode) , 1 r ' ) as f: 
Ig_phi = f . root . IgJgJPoi_phi [ : , 0] 

Jg_phi = f root . IgJgJPoi_phi [ : , 1] 

#Q term indices for phi. For u add +1 to IgQ_phi. 
with open_f ile ( ' IgJgJPQindices_phi%i ' 7, (barcode) r ' ) as f: 
IgQ_phi = f.root IgJgJPQ_phi [ : , 2] 

JgQ_phi = f.root IgJgJPQ_phi [: ,3] 

#Indices for the time derivative Jacobian. 
with open_f ile ( ' IgJgA7»i 1 % (barcode) , 'r' ) as f: 

IgA = f.root IgJgJA[:,0] 
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#Load gradients . 

with open_f ile ( ' Gradients°/ 0 i ' % (barcode) , 1 r ' ) as f: 

Grads = f root . Grads [: ] 

#1 . Compute the G matrix and save it (Most time consuming step). 

boundary _matrix_G(ptot , alphas , Igs ,Ngauss ,pbox, 
normal_vecs , k_save , barcode) 

#Load G matrix. 

with open_f ile("G_matrix_data°/ 0 i" % (barcode) , 'r' ) as f: 

G = f , root , G [ : ] 

Ig_G = array (f . root . IgG [:] ) 

Jg_G = array (f . root . JgG [:] ) 

G=csr_matrix( (G, (Ig_G, Jg_G)) , shape=(5*nq,5*nq)) 

#2. Compute Poisson's stiffness matrix (w/o the 1/2 from IMR) . 

Kg = JPoisson(vol ,ne , Grads, barcode) 

Kg_phi = csr_matrix( (Kg, (Ig_phi , Jg_phi) ) , shape=(5*nq,5*nq) ) 

Kg_u = csr_matrix( (Kg, (Ig_phi+1 , Jg_phi+1) ) , shape=(5*nq, 5*nq) ) 

#2. Introduce the residual s=Gu-I 
Kg_phi [Igs*5+3, Igs*5+3] =-l 
Kg_phi+=G 

#3 . Compute the Q_phi term, (already has the -omega/8 constant) 

Q = JQ_phi (Grads ,ne ,nq,vol , barcode) 

Q_phi = csr_matrix((Q, (IgQ_phi, JgQ_phi)) ,shape=(5*nq,5*nq)) 

Q_u = csr_matrix( (Q , (IgQ_phi+l , JgQ_phi) ) , shape=(5*nq, 5*nq) ) 
Q_BC_u = BC_Q_u(normal_vecs , alphas , areas ,nq) 

#/. Compute the time derivatice term. 

A = fA(alphae,vol,deltat, barcode) 

A = csr_matrix( (A, (IgA, IgA) ) , shape=(5*nq, 5*nq) ) 

#5 . Assemble all constant J terms. Only the Q_u_BC has 
#a minus sign , already included in the function . 

Jfix = (Kg_phi + Kg_u + Q_phi + Q_u + Q_u_BC)*0.5 + A 

#6. Save Jfix. 

data = Jfix. data 

indices = Jfix. indices 

indptr = Jfix.indptr 

with open_f ile ( ' Jf ix_matrix_data_°/ 0 i ' % (barcode) , 1 w ' ) as f: 
atom = Float64Atom(shape=() ) 

Jfix_matrix = f create_vlarray (f root, ' Jfix 1 , atom) 

Jf ix_matrix . append (data) 

Jf ix_matrix . append(indices) 

Jf ix_matrix . append (indptr) 


Algorithm G.l: This function will build the part of the Jacobian J given by the time derivative 
contribution (3.64), the Poisson’s stiffness matrix (5.15) for both the <j> and u and Jacobian with 
respect to m, (3.56). 
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def Jandresidual (x , xO , happ , happO , deltat , normal_vecs , alphas , alphae , 

Igs , alphae , Grads , areas , vol , nq , ne , J , KgphiuG , barcode) : 

II II II 


INPUT: All arguments for all Jacobian and residual functions . 
J is Jfix. 

OUTPUT: A 5nqx5nq csr matrix and a 5nqxl csr column vector. 

II II II 


avg=zeros ( [3 , nq] ) 
dx=x-xO 


#1 . Initiate an array with zeros. 

#2. Compute the difference between configs. 

5] #3. Put each component of mag at 

5] #one row of avg, the natural 

5] #structure for fmalpha. 

dmalpha = fmalpha (avg, alphae) 

x=(x+x0)*0.5 #f. Compute the avg (IMR) config and replace x. 


avg[0] = dx[0 
avg[l] = dx[l 
avg [2] = dx[2 


avg[0] = x[0 

avg [1] = x[l 

avg [2] = x[2 


> : : 5] #5. Build an array where 

. ::5] #each column is the mag 

! : : 5] #at i . 

avgmalpha = fmalpha (avg, alphae) #6. Mag at each (e, alpha). 
avgmalphas = fmalphaS (avg, alphas) #1 . Mag at each (s, alpha). 


avg [0] = happO [0 : 
avg[l] = happO [1: 
avg [2] = happO [2 : 


3] 

3] 

3] 


#8. Do the same for the cte 
#applied field. . . 


avghappalpha = fmalpha(avg, alphae) 

avg[0] = x[3: :5] #9... and the potential phi. 

avgphialpha = fphialpha(avg [0] , alphae) 
avgmtot = fmtot (avgmalpha, ne) 


avg[0] = x0[0 
avg[l] = xO[l 
avg [2] = x0[2 


5] 

5] 

5] 


#10. Malpha structure for mO. 
#No average used. 


malphaO = fmalpha (avg, alphae) 


#Compute all matrix J contributions and residuals . 

J += JD( avghappalpha, avgphialpha, deltat, Grads, 
malphaO , vol , nq , ne , alphae , barcode) 

J += JP_phi (avgmalpha, Grads, nq, vol, barcode) 
r = rA(dmalpha, deltat, nq, vol, barcode) 

r += rD( avghappalpha, Grads, avgphialpha, avgmalpha, dmalpha, 
alphae , vol , nq , deltat , alphae , barcode) 
r += fs_phi_u (KgphiuG, x, avgmalpha, avgmalphas, vol, Grads, Igs, 
nq, alphas ,normal_vecs .areas , barcode) 

B0,B = JEx( avgmalpha, avgmtot, barcode, nq,ne) 

J += B 
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r += rEx (BO ,avgmtot,nq,ne, barcode) 
return J,r 


Algorithm G.2: Assemble all terms in J and r. 
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Appendix H 


Build G matrix 


def boundary _matrix_G(ptot , alphas , Igs ,Ngauss ,pbox, normals ,k_save , barcode) : 

II II II 

Compute all integrals and global indices for G matrix. 

INPUT: vertices locations ,ptot (shape=(3 ,nq) ) ; 

Boundary connectivity matrix, 

alphas, and the indices of boundary vertices 

Igs (column indices of ptot) ; 

Ngauss is the number of vertices within a given triangle s used 
in the integration; pbox=[-xbox, +xbox, -ybox, +ybox, -zbox, +zbox] . 

OUTPUT :3 arrays saved at "G_matrix_data°/i" ’/(barcode) : 

The entries of G (f .root .G) and its 2 corresponding 
global indices for a 5nqx5nq matrix, 

Ig and Jg (f .root . IgG; f .root . JgG) . 

II II II 

"1. Setup up Gaussian locations and weights and the Jacobians" 
epsilon=10** (-15) # Tolerance to detect zero integrals 

xi_eta,weights=gauss_triangle (Ngauss) #Ref triangle Gauss Iocs and weights. 
Jacobians=f J(alphas ,ptot) #Jacobians from the ref triangles . 

Ig, Jg,G= [],[], [] #Where future matrix indices 

# and entries will be kept. 

runonce=False #Create file once. 

print ("Integrals calculations started ") 

"2. Pick a vertex from the boundary and an element 
s and test whether integral is zero" 

for k,i in enumerate (Igs) : #Pick a vertex with 

print("k" ,k, 1 out of ' ,len (Igs) -l)#index i from the surface. 
if k°/„k_save==0 : 

if runonce == False: 

with open_f ile("G_matrix_data7oi" % (barcode) , 'w') as f: 
g_matrix = f . create_earray (f root. 
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1 G 1 , Float64Atom() , shape=(0,)) 

Ig_g_matrix = f . create_earray (f root, 

1 IgG ' , Int32Atom() , shape=(0, ) ) 

Jg_g_matrix = f . create_earray (f root, 

1 JgG ' , Int32Atom() , shape=(0, ) ) 
for g,ig, jg in zip(G, Ig, Jg) : 
g_matr ix . append (g) 

Ig_g_matrix . append (ig) 

Jg_g_matrix . append( jg) 

#Empty the lists. 

G,Ig,Jg= □ , □ , □ 
runonce = True 

print ("File Created and info saved") 

else : 

with open_f ile("G_matrix_data7oi" % (barcode) , 1 a' ) as f: 
g_matrix = f root.G 
Ig_g_matrix = f.root IgG 
Jg_g_matrix = f.root. JgG 
for g,ig, jg in zip(G, Ig, Jg) : 
g_matrix . append ( [g] ) 

Ig_g_matrix . append ( [ig] ) 

Jg_g_matrix . append ( [jg] ) 

#Empty the lists. 

G,Ig,Jg= □ , □ , [] 
print ("saved! ") 

x_i = ptot[:,i] #Location of vertex i. 

for s in range (alphas . shape [1] ) : #Pick a triangular element. 
ig,jg,kg = alphas [:,s] #The 3 vertices locations for s. 

xOs = ptot [ : , ig] 
xls = ptot [ : , jg] 
x2s = ptot [ : ,kg] 

n = normals [:,s] #The normal to element s. 

d=xOs-x_i #The distance from x_i to element s. 

isnormal=dot (n, d) #If isnormal is zero then skip 

if abs (isnormal) < epsilon: # further calculations . 
continue 

"3. Compute the integrals" 

xO,yO,zO = xOs #Build the coordinate transformation 

xl,yl,zl = xls #matrix T from ref triangle 

x2,y2,z2 = x2s #to triangle element s. 

T = array ( [[x2,x0-x2,xl-x2] , 

[y2,y0-y2,yl-y2] , 

[z2,z0-z2,zl-z2]] ) 
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Js = Jacobians [s] #Pick the corresponding Jacobian. 
for a in range(3): #For a given Ns_alpha 

1=0 #compute the integral , I. 

for w,(xi,eta) in zip(weights ,xi_eta) : 
x = dot (T, array ( [1 ,xi , eta] ) ) 
x = x-x_i 

x3=4*pi*norm(x) **3 

f = Ns_a(a,xi , eta) *isnormal* Js/x3 

I+=w*f 

G append ( -I) #(-l) from the derivative 

#of the Green function. 
j = alphas [a] [s] #Append Gif indices fixed by 
Ig. append (i*5+3) #by the Jacobian block structure. 
Jg append (j *5+4) 

"4. Compute the diagonal terms gamma" 
for i in Igs : 

x_i = ptot [ : , i] 

G . append(fgamma(x_i ,pbox) ) 

Ig . append(i*5+3) 

Jg . append(i*5+4) 

with open_f ile("G_matrix_data°/ 0 i" % (barcode) , 'a') as f: 
g_matrix = froot.G 
Ig_g_matrix = f root IgG 
Jg_g_matrix = f root JgG 
for g, ig, jg in zip(G,Ig, Jg) : 
g_matrix . append ( [g] ) 

Ig_g_matr ix . append ( [ig] ) 

Jg_g_matrix . append ( [jg] ) 
print ("Last save!") 


Algorithm H.l: Algorithm to compute G. 
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Appendix I 


Update Functions 


def Update (xO , hapO , k_lim , k_save , args , args_phi , 

args_u,args_gmres .barcode ,barcode_hist) : 

II II II 

INPUT: xO is the initial configuration with shape=(5*nq, ) ; 

hapO is the magnetic field with shape=(3*nq, ) ; k_lim is the number 

of time steps; k_save gives the steps to save xO; 

args is a dictionary with info for Newton_Raphson function; 

args_phi , args_u, args_gmres are the arguments 

for the GMRES algorithm precontioned with AMG for 

phi and u blocks; barcode is the name of the file where we 

saved the mesh info (ex: indices for the residual and J); 

barcode_hist is the name of the file where we save x. 

OUTPUT: a file "xvst°/,i" %(barcode_hist) with an array 
with shape=(timesteps ,5*nq) with timesteps=k_lim%k_save. 

II II II 

save_j=0 # Current save number=row of "xvst" file. 

runonce=False #To create file once. 

#l.For a given time step k. 
for k in range (k_lim) : 

#2. Test whether current k is a multiple of k_step to save. 
if k°/„k_save==0 : 

if runonce==False : 

with open_f ile("xvst%i" °/ 0 (barcode_hist) , 1 w 1 ) as f: 
steps=k_lim/k_save 
shape_x= (steps , 5*nq) 
atom=Float64Atom() 

history_x=f create_carray (f root x 1 , atom, shape_x) 
history_x [save_j] =x0 
save_j+=l 
runonce = True 

else : 
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with open_f ile("xvst%i" °/„(barcode_hist) , 'a') as f: 
f root . x [save_j] =x0 
save_j+=l 

#3 . Initiate time counter to predict cumputational time 
#for the k_lim iterations. 
tO=perf _counter () 

#f.Set the first particular configurations as xO. 
x=array (xO) 

xO=Newton_Raphson (x , xO , hapO , hapO , args , args_phi , 
args_u , args_gmres , barcode) 

tl=perf _ count er () 

#5. Print time left until k_lim. 

print("NR #°/„i started out of '/,i, took °/ 0 .4f secs" 7o(k,k_lim, (tl-tO) ) ) 

print (" time left: °/ 0 .4f hrs" °/.( (k_lim-k) * (tl-tO) /3600) ) 

k += 1 

print ("UPDATE FINISHED ! ! ! ") 


Algorithm 1.1: The function that computes the evolution of the configuration of the system from initial 
conditions. 

def Newton_Raphson (x , xO , hap , hapO , args , args_phi , args_u , args_gmres , barcode) : 

II II II 

INPUT: Arguments decribed in Update function. 

OUTPUT: The congiguration x associated to t+deltat . 

II II II 

#O.Get variables out of args. 

deltat=args . get ( 1 deltat ' ) 
normal_vecs=args . get ( ' normal_vecs 1 ) 
alphas=args . get ( 1 alphas ' ) 
alphae=args . get ( 1 alphae ' ) 

Igs=args . get ( ' Igs 1 ) 
alphac=args . get ( 1 alphae ' ) 

Grads=args . get ( ' Grads ' ) 
vol=args . get ( ' vol 1 ) 

Jf ix=args . get ( ' Jf ix 1 ) 

KgphiuG=args . get ( 1 KgphiuG 1 ) 
it_max=args . get ( 1 it_max ' ) 
restol=args . get ( 1 restol ' ) 

#1 . Initiate sequence of linear systems. 
it=0 

while 1 : 

#2. Build Jacobian matrix and the residual . 

J,r = Jandresidual(x,xO, hap, hapO, deltat, normal_vecs, 

alphas , alphae , Igs , alphae , Grads , areas , 
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vol ,nq,ne , Jfix,KgphiuG, barcode) 

#3 .Compute residual norm. 

normres=norm(r) 

print (" NR it=" , it , "norm res" ,normres) 

#f.Test whether residual is small enough or no its is exceded. 
if it>it_max or normres<restol : 
return x 

#5. Solve a linear system of equations . 

x+=rILU_AMG_GMRES (J, -array (r) reshape (5*nq, ) , 

array (xO) , args_phi , args_u, args_gmres) 

it+=l 


Algorithm 1.2: Our version of the Newton-Raphson algorithm with GMRES with preconditioners ILU 
and AMG (rILU_AMG_GMRES is not described in this thesis). 
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Appendix J 


Amplitude and Phase Distribution 
for Standard Problem 


From the OOMMF magnetization evolution data available from [27] we have the following 
amplitude and phases distributions for m x , m y and m z for the two resonance peaks: 8.25 
GHz and 11.25 GHz. 
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2.25 



Figure J.2: OOMMF and 11.25 GHz. 
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Appendix K 


Discrete 


Fourier Transform 


We start with the Fourier series and obtain the approximation to the integrals of the coeffi- 
cients using the trapezoid rule following [12]. Then given them a matrix representation. The 
goal is to understand what this operation so that the DFT data computed in Chapter 6 is 
understood. 

Consider a periodic function or a slice of a function that is then copied along the x axis. 
The representation of this function is given by 


/(*) 


E 


i2Tckx 

Cfce l 


k =— oo 




(K.l) 

(K.2) 


In practice the function might be too complicated to integrate or discrete if it results from 
a laboratory experiment or numerical simulations. This implies that the integrals have to be 
approximated numerically. For that we can use the Trapezoid Rule [12]. This method has 
the purpose of approximating the integral as follows. First we divide the a generic interval 
of integration [a, b] into N segments of index k. The area above each segment is given by 
the average of the values of f{x) at each point that bounds the segment k times the width of 
the segment, h = L/N. The integrals is then approximately 

r b N 

/ f(x)dx « ^2 A e (K.3) 

Ja 6=1 


where the higher is the number N the closer is the A e to the actual area under f(x) 

A e = ,, /(“+ (e-lW +/(“ + <*> (K.4) 

The point on the left of segment e is Xq = a+(e— l)h while the point on the right is xf = a+eh. 

Observe that the area above each segment is computed with information associated with 
each point that determined it. The integral is then given substituting (K.4) into (K.3), 
isolating the first and last points of [a, b\ the integral is 
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f f{x)dx k, h 
J a 


N—l 


1 /(«) + \f(b) + ^ /(a + eh) 

e=l 


(K.5) 


2 2 " 

the sum goes through all points on the right of each segment and that are not on the boundary. 

We can use now the result (K.5) to approximate the calculation of the coefficients, we 
obtain then 
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(K.6) 


where n is the index of the point within the ]o, b[ and x n = 0 + nh. Observe that since f(x) 
is periodic then /( 0) = f(L ) hence 
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(K.7) 


n= 0 


The first term of the sum n = 0 results from the average between the f(x n ) at the two 
boundary points x n = a and x n = b. 

The calculation is done summing N areas associated with N segments that divide [0,L], 
Each area depends on two points that define the segment. Since the function is periodic the 
value of the function at the boundary points is equal hence we lump them leaving N — l 
calculation in the sum in (K.5). This average is then introduced in the sum as n = 0. 

The transformation (K.7) can be used starting with a continuous function which is eval- 
uated at N points within [xo,xn~i] where xn = xq were left or we could start with data 
already discretized. 

With the N coefficients c& how do we invert the expression and obtain the values of f(xk) 
again? Multiplying both sides of (K.7) by e n and summing for all k we have 
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(K.8) 


where we substituted the definition of h and renamed f n = f(x n ). On the RHS we observe 
the second sum is a geometric series that is equal to 


N-l 

^ ^ gi27rkm/N g—i2nkn/N 

k = 0 


Subtituting in (K.8) we have finally 




N-l 

fm — ^ " c^.e 

k= 0 


(K.9) 


(K.10) 


Observe from the N coefficients we can recover all original discrete f(x m ) values, there is 
no loss of information. 

If we define the fundamental constant oj = e l2n / N , we can organize the mapping of coeffi- 
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cients in (K.10) into all f m values as the following matrix operation 
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(K.ll) 


Observe that (K.10) gives how a particular entry of f on the LHS of (K.ll) is equal to a dot 
product of a row of the matrix Fy and with the vector c. Notice also the matrix is symmetric. 
To obtain the vector of coefficients, c, from the values f we have 


c — 


(K.12) 


The key to compute such an inverse is to know the relation among the basis vectors in the 
columns of Fy. If we choose two columns n and mn of Fy and compute the dot product 
observe what we are actually computing is the expression (K.9). When n = m we get N, the 
squared norm of the vector, if n ^ m then we have 0, meaning if we normalize all vectors 
dividing yFy we have a orthonormal basis 
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Vn 


F; 
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Fy = / 


(K.13) 


where we transpose and take the complex conjugate of the matrix. This matrix relates the 
basis vectors just as (K.9) and from it we can obtain the inverse yielding 


Fn = yh? (K.14) 

Hence expression (K.12) is given by 

c = ^Fyf (K.15) 

where we used the fact that the matrix is symmetric, Fy = Fy. Notice that entry k in c is 
given by a dot product of row k with the f vector, that is the expression (K.7). The meaning 
of such an operation can better be observed from an example. Suppose f corresponds to data 
for N time steps and the variable x is replaced by time, t. Then expression (K.7) can be 
written as follows 

N—l 

c k^^T,f^ e ~ i2nfktn ( K - 16 ) 

71=0 

where we introduce the following definitions 


fk 


k 

NAt 

nAt 


(K.17) 

(K.18) 


The sequence of entries in a row corresponds to evaluating e * 27r F*« at a sequence of instants 
of time starting at to = 0 up to fy_i = (N — 1)A t. Its real part or imaginary part represent a 
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discretized wave with frequency /*,. The bottom rows of F]y (or Fn) have higher frequencies 
that the top ones. Now the vector f is decomposed in a linear combination of these vectors 
with coefficients which are complex numbers obtained by computing the dot product, that 
is by projecting f into each particular oscillation. The coefficient c * is characterized by two 
numbers, the amplitude and phase, qualitatively with both we measure how much f ’’looks 
like” each of these basis vectors. 
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Appendix L 


Software Used 


• The code is written using: Anaconda Software Distribution. Computer software. Ver. 
4.3.1. Continuum Analytics, Nov. 2017. Web. https://continuum.io. 

• The pictures in Chapter 3 and 4 were made with: Harrington, B. et al (2004-2005). 
Inkscape. http:/ /www. inkscape.org/. 
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